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Abstract 

Current  tracking  and  adaptive  optics  techniques  cannot  compensate  for  fast- 
moving  extended  objects,  which  is  important  for  ground-based  telescopes  providing 
space  situational  awareness.  To  fill  this  need,  a  vector-projection  maximum-likelihood 
wave-front  sensing  algorithm  development  and  testing  follows  for  this  application.  A 
derivation  and  simplification  of  the  Cramer-Rao  Lower  Bound  for  wave-front  sensing 
using  a  laser  guide  star  bounds  the  performance  of  these  systems  and  guides 
implementation  of  a  vastly  optimized  maximum-likelihood  search  algorithm.  A  complete 
analysis  of  the  bias,  mean  square  error,  and  variance  of  the  algorithm  demonstrates 
exceptional  performance  of  the  new  sensor.  A  proof  of  concept  implementation  shows 
feasibility  of  deployment  in  modern  adaptive  optics  systems.  The  vector-projection 
maximum-likelihood  sensor  satisfies  the  need  for  tracking  and  wave-front  sensing  of 
extended  objects  using  current  adaptive  optics  hardware  designs. 
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One  More  Roll 

We  toast  our  faithful  comrades  now  fallen  from  the  sky 
And  gently  caught  by  God’s  own  hand  to  be  with  him  on  high. 

To  dwell  among  the  soaring  clouds  they  knew  so  well  before 
From  dawn  patrol  and  victory  roll  at  heaven ’s  very  door. 

And  as  we  fly  among  them  there  we  ’re  sure  to  hear  their  plea 
“Take  care,  my  friend,  watch  your  six,  and  do  one  more  roll...  just  for  me.  ” 

Gerald  (Jerry)  Coffee,  Captain,  USN  (Ret.) 
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MULTI-DIMENSIONAL  WAVE  FRONT  SENSING  ALGORITHMS  FOR 
EMBEDDED  TRACKING  AND  ADAPTIVE  OPTICS  APPLICATIONS 


I.  Introduction 

Atmospheric  turbulence  affects  clarity  of  anything  in  space  viewed  through  large 
telescopes.  Machines  that  perfonn  optical  tracking  of  moving  targets  or  provide  high- 
resolution  imaging  must  correct  for  turbulence  effects  by  detecting  the  distorted  wave- 
front  caused  by  turbulence  to  prevent  loss  of  tracking  ability  or  image  corruption  [16]. 
The  capability  to  detect  distortion  in  the  wave-front,  or  relative  position  changes  of  an 
image,  is  often  embedded  or  built  into  the  wave-front  sensor  and  processing  algorithms 
as  a  part  of  the  adaptive  optics  system  [16].  Modern  adaptive  optics  systems  allow  for 
wave-front  correction  with  only  guide  stars  and  small,  extended  sources,  sometimes 
requiring  post-processing  of  gathered  data  [12].  A  new  maximum-likelihood  wave-front 
sensing  algorithm  embedded  in  proven  adaptive  optics  designs  could  enhance  detection 
for  non-ideal  conditions  and  real-time  operations  [5].  What  and  where  atmospheric 
turbulence  is,  how  adaptive  optics  attempt  to  overcome  the  effects  of  this  turbulence,  and 
why  these  optics  systems  need  improvement  all  become  clear  in  the  following  sections. 
1.1.  Background  and  Motivation 

1.1.1.  The  Effects  of  Atmospheric  Turbulence 

Many  factors  on  earth,  such  as  natural  processes  and  terrain  features,  affect 
weather  significantly;  however,  the  main  driver  of  turbulence  is  the  sun’s  uneven  heating 
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of  the  earth’s  surface.  The  uneven  heating  causes  convection  currents  and  wind 
spawning  circular  currents,  eddies,  which  trap  varying  temperatures  throughout  the 
atmosphere  causing  variations  of  the  index  of  refraction  thereby  distorting  the  wave- 
front.  Figure  1  shows  the  first  two  major  layers,  the  troposphere  and  stratosphere, 
containing  99.9  %  of  the  earth’s  atmosphere,  and  whose  turbulence  is  responsible  for  the 
majority  of  light  distortions  [13].  The  figure  also  indicates  an  average  temperature 
gradient;  a  few  realistic  sample  temperature  gradients  as  seen  through  different  columns 
of  air;  and  other  sources  of  turbulence  such  as  shearing  winds,  terrain,  and  natural 
processes  feeding  convection.  The  results  of  these  sources  of  turbulence  can  combine  to 
distort  a  wave-front  as  it  passes  through  different  temperature  gradients  in  the  earth’s 
atmosphere. 
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Figure  1.  Temperature  Gradients  and  Turbulence  Sources  in  the  Atmosphere  [13] 
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A  clearer  view  of  how  a  wave-front  distorts  and  how  the  wave- 
front  initially  forms  is  available  in  Figure  2.  A  point  source,  or  a 
distant  star,  emits  light,  which  travels  outward  from  the  star  much  as 
ripples  travel  outward  from  a  pebble  thrown  in  a  pond.  When  these 
“waves”  are  far  away  from  the  source,  they  appear  as  a  straight  line, 
forming  a  wave-front.  Researchers  often  model  the  propagating  waves 
as  a  two-dimensional  Fourier  Transform.  Much  like  taking  the  Fourier 
Transform  of  a  single  point  in  time  results  in  a  straight  line  in  the 
frequency  domain,  a  point  source  in  space  transforms  to  a  plane  wave 
related  to  spatial  frequency  rather  than  temporal  frequency  [8].  The 
wave-front  does  not  distort  much  as  it  passes  through  the  stratosphere, 
as  temperature  variations  seldom  occur  there;  however,  the  troposphere 
severely  distorts  the  wave-front  due  to  the  numerous  opportunities  for 
eddies  to  form  and  trap  temperature  variations.  The  result  is  a 
corrupted  wave-front  that,  when  focused  onto  an  imaging  device 
produces  a  blurry  and  distorted  image  bearing  little  resemblance  to  the 
original  object. 

In  addition  to  using  phase  screens  to  model  the  temperature 
variations  and  the  relative  refractive  index  changes  directly  at  different 
altitudes,  researchers  use  Zernike  polynomials,  or  “Zemikes”,  to 
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Figure  2. 
Distorting  of 
Wave-Front 
Moving 
Through 
Atmosphere 
[12] 


characterize  the  distortions  in  the  wave-front  itself  [16].  As  opposed  to  a  rectangular 
based  set  of  polynomials,  the  Zernike  polynomials  describe  a  set  of  circular-based,  two- 
dimensional  functions  corresponding  to  the  circular  opening  in  a  telescope  or  other 
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imaging  device  [16,  22].  These  models  are  crucial  to  correcting  the  wave-front  in  an 
Adaptive  Optics  (AO)  system. 

1.1.2.  Adaptive  Optics  Solutions 

Although  there  are  many  applications  for  adaptive  optics  in  modern  imaging 
systems,  the  basic  structure  as  shown  in  Figure  3  for  a  general  large  telescope  system 
remains  relatively  constant  across  the  applications  [16]. 

Light  From 
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A  simple  trace  through  the  system  reveals  that  light  enters  through  the  telescope  lens  with 
a  distorted  wave-front,  and  then  reflects  from  an  adaptive  mirror,  a  mirror  that  can 
deform  using  mechanical  actuators,  which  is  initially  flat  as  there  is  no  information  to 
correct  the  wave-front.  The  light  then  continues  to  a  fifty-fifty  beam-splitter  sending  half 
of  the  light  into  a  lens,  which  focuses  the  light  onto  a  high-resolution  imaging  device,  and 
the  other  half  to  the  wave-front  sensor.  The  first  portion  of  the  wave-front  sensor  both 
optically  and  electrically  detects  measurable  parameters  of  the  wave-front,  passing  that 
information  to  an  algorithmic  portion  of  the  wave-front  sensor  to  estimate  the  parameters 
for  later  modeling.  Since  the  wave-front  sensor  is  the  heart  of  this  system,  this  thesis 
concerns  itself  with  the  algorithmic  portion  of  the  wave-front  sensor.  These  estimates 
pass  to  the  reconstructor  in  the  control  system,  which  builds  a  model  of  the  wave-front 
and  then  applies  that  information  to  a  known  model  for  the  adaptive  mirror  to  attempt 
wave-front  correction.  If  the  wave-front  has  a  lag  or  dip  in  it  causing  the  light  to  arrive 
later  than  expected,  the  mirror  must  have  a  corresponding  bump  to  accelerate  the  light 
back  to  its  appropriate  phase  to  compensate  for  the  distortion.  Once  an  initial  estimate 
corrects  the  wave-front,  additional  wave-front  sensing  refines  the  current  estimates  and 
detects  further  changes,  producing  higher  quality  results  for  future  images. 

1.1.3.  Problems  and  Need  for  Improvement  in  Wave-Front  Sensing 

The  applications  for  higher  quality  imaging  span  the  gamut,  from  ground-based 
and  space-based  telescopes  to  military  applications  such  as  the  Airborne  Laser  and  even 
medical  services  such  as  measuring  aberrations,  or  deformities,  in  an  eye.  These 
applications  drive  the  need  for  better  quality  imaging  and  improvements  in  wave-front 
sensing. 
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As  with  any  scientific  research,  improvements  require  a  metric  by  which  to 
measure  results  and  draw  conclusions.  To  this  end,  knowing  the  structure  for  a  Cramer- 
Rao  lower  bound  (CRLB)  would  provide,  independent  of  the  estimation  technique,  an 
analytical  method  to  judge  the  efficacy  of  current  and  proposed  wave-front  sensing 
algorithms.  Once  known,  the  Cramer-Rao  lower  bound  can  also  guide  research  for 
improving  current  estimation  techniques  as  well  as  developing  new  estimation 
approaches  to  manage  more  complex  imaging  scenarios. 

A  complex  situation  of  interest  is  imaging  of  extended  objects,  or  light  sources 
that  do  not  conform  to  the  definition  of  a  point  source,  such  as  a  satellite  in  orbit,  the 
surface  roughness  of  the  sun,  a  scud  missile,  or  even  a  truck  on  a  highway.  Tracking  a 
satellite  in  orbit  allows  for  space  situational  awareness,  or  imaging  of  foreign  assets, 
without  placing  costly  assets  in  space;  however,  it  requires  wave-front  updates  for  this 
extended  object,  the  satellite,  at  an  incredible  rate  of  1,000  Hz  or  greater  due  to  the  speed 
in  which  the  satellite  moves.  A  complication  to  the  satellite -tracking  scenario  stems  from 
the  typical  optical  tracking  system,  which  causes  the  image  to  fill  the  field  of  view  and 
allows  new  information  to  enter  the  scene  while  tracking,  defeating  current  fast-acting 
sensors.  The  surface  intensity  variation,  or  roughness,  of  the  sun  is  a  unique  problem  in 
that  the  image  gathered  has  extremely  low-contrast  features  for  tracking  or  correlation, 
disabling  most  modern  wave-front  sensors;  but  imaging  of  the  sun  is  necessary  to  predict 
communication  outages  and  solar  weather  in  general.  The  scud  missile,  truck,  and  other 
daytime  AO  applications  represent  a  class  of  objects  whose  backgrounds,  like  the  sun, 
are  not  black  reducing  the  contrast,  and  require  rapid  wave-front  updates  for  the  dynamic 
turbulence  between  the  object  and  imaging  device.  Imaging  extended  objects,  although 
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merely  a  collection  of  point  sources  with  spatial  reference  to  each  other,  is  not  the  only 
type  of  imaging  that  current  wave-front  sensors  can  have  poor  perfonnance. 

Occasionally,  atmospheric  turbulence  results  in  a  tip-tilt,  represented  by  Zernike 
polynomials  two  and  three,  beyond  the  physically  measurable  range  of  the  sensor  causing 
an  unknown  in  the  collection  of  estimated  parameters  and  preventing  the  reconstructor 
from  modeling  the  wave-front.  This  unknown  occurs  when  the  tip  or  tilt  is  so  great  that 
the  majority  of  the  image  moves  off  the  detector  leaving  the  algorithm  a  small  amount  of 
information  to  work  with.  Modern  sensors  are  not  capable  of  controlling  such  a  situation, 
and  the  entire  adaptive  optics  system  suffers  when  a  single  sensor  cannot  acquire  an 
accurate  estimate  for  the  wave-front. 

Although  the  optical  and  electrical  properties  of  current  sensors  potentially 
support  the  previously  mentioned  improvements,  the  algorithms  currently  in  use  do  not; 
therefore,  an  investigation  of  a  vector-projection,  maximum-likelihood-correlating  wave- 
front  sensor  guided  by  Cramer-Rao  lower  bounds  and  simulation  experiments  will 
proceed. 

1.2.  Summary  of  Current  Techniques 

Several  factors  limit  the  performance  of  current  adaptive  optics  techniques 
preventing  the  ability  to  track  or  perform  wave-front  sensing  for  fast-moving,  extended 
objects,  or  low-contrast  objects.  The  largest  contributor  to  these  limitations  is  the  wave- 
front  sensor,  which  provides  the  necessary  information  for  the  adaptive  optics  system  to 
correct  the  wave-front  deformities. 

There  are  two  main  categories  of  wave-front  sensors,  low-order  wave-front 
sensors  and  advanced  wave-front  sensors,  both  of  which  are  capable  of  detecting  wave- 
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front  distortions.  The  advanced  wave-front  sensors  typically  produce  better  performance 
through  higher  order  computations  and  more  complex  algorithms;  however,  most  cannot 
image  an  extended  object  and  none  are  capable  of  the  tracking  application  as  closed  loop 
speeds  are  currently  very  low  [6].  Current  low-order  wave-front  sensors  provide  slightly 
lower  imaging  performance,  but  operate  at  up  to  1000  Hz,  allowing  for  tracking  and  other 
fast-moving  imaging  applications  [16].  These  simpler  estimation  techniques  include 
numerous  wave-front  sensors;  however,  only  the  easily  implemented  and  fast-operating 
Shack-Hartmann  and  Short- Wavelength  Adaptive  Techniques  (SWAT)  wave-front 
sensors  are  common  today  [16].  Both  of  these  sensors  use  a  centroid-based  algorithm  to 
estimate  tip  and  tilt,  and  this  algorithm  can  have  extremely  poor  performance  when 
attempting  wave-front  sensing  or  tracking  on  an  extended  object  [16].  The  simplicity  of 
the  centroid  algorithm  suggests  that  a  more  complex  and  statistically  based  algorithm 
could  surpass  these  sensors  in  performance,  possibly  retaining  the  operating  speed  while 
tracking  or  performing  wave-front  sensing  on  extended  objects. 

Research  indicates  the  theoretical  vector-projection  maximum-likelihood  wave- 
front  sensor  can  achieve  the  performance  of  a  low-order  wave-front  sensor  for  tracking 
and  wave-front  sensing  of  guide  stars  while  providing  suitable  performance  for  imaging 
extended  objects  [5].  This  wave-front  sensor  uses  the  same  hardware  system  as  the 
SWAT  wave-front  sensor;  however,  the  algorithm  is  a  maximum-likelihood  estimation 
technique,  which  provides  correlation  capability  for  an  extended  object  while 
maintaining  perfonnance  for  point  sources  [5].  Currently  only  limited  simulated 
statistical  characterization  of  this  sensor  is  available  and  the  tracking  application  requires 
a  modern  processor  to  implement  this  more  complex  algorithm  [5]. 
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1.3.  Contributions  and  Scope 

It  is  the  goal  of  this  research  to  quantify  the  efficacy  of  a  vector-projection, 
maximum-likelihood-correlating  wave-front  sensor  for  tracking  extended  objects  based 
on  a  satellite  application,  as  well  as  a  couple  wave-front  sensing  applications,  through 
three  facets  [5]. 

The  first  contribution  is  generalized  model  for  the  Cramer-Rao  lower  bound  with 
assumptions  allowing  for  future  applications  provides  the  analytical  basis  for  research. 
The  CRLB  should  be  applicable  to  any  type  of  wave-front  sensor. 

The  second  contribution  is  an  algorithmic  analysis  to  increase  the  temporal 
performance  of  the  new  complex  maximum-likelihood  algorithm  to  allow  simulations 
that  thoroughly  characterize  the  noise-independent  bias  of  the  algorithm  resulting  in  a 
third  contribution  as  well  as  the  noise  statistics  of  mean  squared  error  (MSE)  and 
variance  (VAR)  for  a  fourth  contribution.  The  variance  directly  compares  to  the  Cramer- 
Rao  lower  bound  to  reveal  limitations  in  the  algorithm.  The  search  algorithm  developed 
in  this  phase  contributes  to  other  applications  requiring  a  fast  and  complete  algorithm  to 
perform  the  search  of  functions  with  special  properties  such  as  maximum-likelihood. 

The  fifth  contribution  is  a  proof  of  concept  for  a  feasible  method  of  developing 
this  algorithm  for  embedded  hardware  implementation  and  a  complete  plan  for 
implementation  offers  insight  to  researchers  in  the  field  looking  for  feasible  solutions. 
The  third  criterion  is  complete  when  a  single  working  implementation  emerges;  however, 
multiple  revisions  provide  further  utility. 
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This  three-faceted  exploration  secures  a  concrete  approach  to  the  research, 
development,  and  implementation  of  a  vector-projection,  maximum-likelihood- 
correlating  wave-front  sensor. 

1.4.  Approach/Methodology 

The  three-faceted  investigation  above,  with  the  provided  motivation,  is  a  template 
that  guides  the  organization  of  both  the  research  and  this  document.  A  thorough 
investigation  of  current  techniques  with  appropriate  discussion  of  relevant  subjects 
provides  the  necessary  foundation  for  research.  This  leads  to  development  of  the 
tracking  and  wave-front  sensing  application  environments  for  producing  realistic 
simulations  and  allowing  accurate  characterization  of  the  new  algorithm.  Theoretical 
analysis  develops  the  CRLB  for  wave-front  tilt  estimates,  which  provides  input  for 
development  of  a  fast,  compact,  and  complete  search  algorithm  for  discovering  the  peak 
likelihood.  From  the  validated  algorithm  extends  a  focused  hardware  implementation. 
The  results  of  extensive  simulations  provide  the  bias,  mean  squared  error,  and  variance 
statistics  characterizing  the  algorithm  for  tracking  and  numerous  wave-front  sensing 
applications.  The  research  concludes  with  a  synopsis  and  areas  of  further  research, 
allowing  for  future  contributions  to  the  body  of  knowledge  regarding  tracking  and  wave- 
front  sensing. 
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II.  Background 


2.1.  Current  Wave-Front  Sensing  Limitations 

Modern  imaging  of  extended  objects  requires  either  a  stable  point  source  in  the 
field  of  view  or  complex  optics  and  algorithms  to  detect  the  wave-front  correctly  across 
the  lens  of  the  telescope.  The  extended  source  typically  forces  researchers,  astronomers, 
and  field  users  to  find  or  create  a  guide  star  close  to  the  extended  object  they  wish  to 
view.  The  applications  mentioned  previously,  particularly  imaging  large  objects  or 
tracking  fast  moving  targets,  are  difficult  or  impossible  to  realize  using  nearby  or 
artificial  guide  stars.  Wave-front  sensing  and  tracking  is  possible  due  to  the  complex 
system  mentioned  in  the  introduction;  however,  the  key  components  are  the  wave-front 
sensor  and  the  algorithm  to  detennine  tip  and  tilt.  The  following  describes  the  typical 
tip-tilt  only  detectors  and  a  few  more  complex  wave-front  detection  methods,  directly 
compares  and  summarizes  the  features  of  each  sensor,  and  finally  presents  areas  of 
potential  research  given  this  information. 

2.2.  Low  Order  Wave-Front  Sensors 

2.2.1.  Shack-Hartmann  Wave-Front  Sensor 

The  most  widely  employed  wave-front  sensor  uses  the  Hartmann  test  to  estimate 
the  linear,  lower  order  Zernikes,  two  and  three,  and  currently  has  the  best  overall  real¬ 
time  performance  [16].  Both  the  hardware  structure  and  the  algorithm  to  gather  offset,  or 
tip  and  tilt,  information  lead  to  a  simple  mathematical  model  stemming  from  the 
elementary  nature  of  the  sensor  as  described  below  [17]. 
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Figure  4  illustrates  a  typical  Shack-Hartmann  sub-aperture  array  in  one- 
dimension  and  a  single  sub-aperture  in  Figure  5  indicates  a  linearly  tilted  wave-front  and 
the  corresponding  offset  in  two-dimensions  when  focused  [16].  The  sub-apertures  must 
be  small  enough  to  meet  the  Nyquist  sampling  criterion  to  ensure  that  the  curved  wave- 
front  is  linear  in  the  region  measured  by  the  sensor  driving  the  overall  number  of  sub¬ 
apertures  [8],  The  Nyquist  rule  applies  to  any  sub-aperture  type  system,  as  well  as 
another  generalized  rule  that  imposes  a  requirement  of  approximately  one  adaptive  optics 
channel,  sub-aperture,  per  turbulence  coherence  radius  r0  as  a  minimum,  independent  of 
the  telescope  size  [14],  Larger  numbers  of  sub-apertures  implies  smaller  sizes;  however, 
this  larger  number  of  sensors  can  introduce  more  noise  into  the  system  and  decrease 
light,  degrading  overall  system  perfonnance  [5]. 


Incident 
Optical  Field 


Figure  4.  Wave-Front  Sensor  Array  [16]  Figure  5.  Single  Wave-Front  Sensor 

Element  [16] 
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These  sub-apertures  consist  of  a  lenslet  array,  which  focuses  the  light  onto  a  charge 
coupled  device  (CCD)  array  for  the  individual  wave-front  sensors  (WFS)  [16].  The  CCD 
readout,  where  the  information  collects,  is  the  second  opportunity  for  significant  noise 
injection  before  the  wave-front  algorithm  begins  processing. 

The  algorithm  driving  a  Shack-Hartmann  wave-front  sensor  is  a  simple  two- 
dimensional  centroiding  algorithm  [17].  Each  intensity  readout  multiplies  a  linear 
position  number,  then  average  together,  and  finally  the  total  power  in  the  image  divides 
the  result  for  the  centroid  in  one-dimension  and  then  repeats  for  the  next  dimension.  This 
operation  takes  a  minimal  amount  of  time  allowing  greater  than  100  Hz  operation,  and 
provides  quality  results  for  guide  stars  and  moderately  extended  objects  [12].  There  is  a 
lower  bound  on  error  for  shot  noise,  or  quantization  noise;  however,  it  is  somewhat 
restrictive  and  only  applies  to  the  Shack-Hartmann  wave-front  sensor  and  guide  stars 
[16].  The  simplistic  nature  of  this  algorithm  lends  itself  to  improvement  in  accuracy  as 
time  permits  such  investigations. 

One  performance  improvement  for  the  Shack-Hartmann  sensor  came  from 
research  at  MIT  Lincoln  Laboratory;  the  short  wavelength  adaptive  techniques  wave- 
front  sensor,  which  splits  the  incoming  light  to  two  lenslet  arrays  and  two  CCDs  oriented 
at  90°  to  each  other.  The  performance  improvement  stems  from  allowing  the  CCD  to 
gather  all  of  the  charge  in  the  image  into  vector  readout,  or  a  projection,  and  then 
performing  a  one-dimensional  centroiding  algorithm  for  each  orientation  [2].  Although 
image  projection  allows  for  both  faster  readout  and  lower  readout  noise,  it  decreases  the 
brightness  of  the  original  image,  decreasing  the  signal  to  noise  ratio  (SNR)  making  an 
accurate  estimate  less  likely. 
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2.2.2.  Shearing  Interferometer 

A  more  complex  wave-front  sensor  not  typically  considered  outside  of  academia 
that  strictly  estimates  Zemikes  two  and  three  is  the  lateral  shearing  interferometer  [23], 
Although  the  physical  implementation  of  a  single  shearing  interferometer  can  be  simple, 
the  algorithm  to  retrieve  a  usable  tip-tilt  requires  a  high  degree  of  effort,  and  the  model 
for  the  wave-front  sensor  system  clearly  indicates  the  non-mathematical  foundation  of  the 
apparatus  and  the  amount  of  processing  required  to  retrieve  phase  information  [16]. 

The  physical  apparatus  splits  the  incoming  light  several  times  encompassing  the 
entire  wave-front  of  the  sensor  to  perfonn  filtering  and  polarization  for  different 
measurement  techniques  [16].  Once  split,  the  beam  splits  again  before  shearing  in 
orthogonal  directions  by  a  tunable  amount,  only  to  recombine  with  the  non-sheared 
version  and  create  an  interference  pattern  as  shown  in  Figure  6  [23],  This  pattern  has  a 
direction  relationship  to  the  wave-front  tilt,  and  a  sinusoidal  nature  over  time  allowing 
researchers  to  correct  the  wave-front  in  a  reasonable  time  frame  [16]. 


Distorted 

Wavefront  Y  shear 


Figure  6.  Shearing  Interferometer  Final  Stage  Operation  [16] 
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Decoding  the  phase  from  this  interference  pattern  takes  many  forms;  however,  all 
algorithms  lead  to  similar  results  with  a  modest  time  delay  and  correct  operation  for  point 
sources  [16].  The  limitation  to  point  sources  stems  from  the  expectation  of  a  plane  wave 
at  the  receiver.  Without  a  point  source,  the  interference  pattern  includes  noise  from  the 
shape  of  the  extended  object  and  corrupts  the  output  waveform.  The  lateral  shearing 
interferometer  is  the  most  tunable  wave-front  sensor,  but  tuning  is  crucial  to  match  the 
Shack-Hartmann  sensor  under  ideal  conditions. 

2.3.  Advanced  Wave-Front  Sensors 

2.3.1.  Curvature  Wave-Front  Sensor 

A  promising  new  wave-front  sensor  is  the  curvature  wave-front  sensor. 

Curvature  sensing  has  additional  requirements  for  the  adaptive  optics  system  by  adding  a 
secondary  deformable  mirror  [1].  The  hardware  for  this  system  relies  on  the  Shack- 
Hartmann  or  other  low-order  wave-front  sensing  detectors;  however,  the  algorithm 
driving  the  higher  order  results,  Zernikes  four  and  above,  takes  the  same  information  and 
performs  a  superior  analysis  at  an  elevated  processing  cost  [1].  The  key  for  this  method 
is  the  requirement  for  an  accurate  tip-tilt  sensor  in  order  to  perform  correctly,  thus 
requiring  the  best  low-order  Zemike  sensor/estimator  possible. 

Aside  from  the  addition  of  a  defonnable  mirror  shown  in  Figure  7,  the  first  mirror 
corrects  for  tip-tilt  and  the  second  correct  higher  order  Zernikes,  the  fundamental  concept 
of  sampling  the  image  changes  for  a  curvature  wave-front  sensor  [1].  Sampling  of  in 
focus  and  out  of  focus  images  occurs  simultaneously  at  a  minimum  of  1  kHz  rate  using  a 
special  parabolic  mirror  as  in  Figure  8. 
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Figure  7.  Curvature  Adaptive  Optics  Setup  Figure  8.  Curvature  Sensing  Setup  [1] 
[1] 


The  multi-phase  sampling  allows  wave-front  correction  over  the  entire  visible  spectrum, 
and  provides  the  flexibility  to  operate  at  lower  frequencies  as  well;  however,  like  the 
shearing  interferometer  it  assumes  a  point  source  is  the  subject  of  the  image  [1].  The 
complexity  of  this  system  forces  the  researcher  to  justify  the  modest  performance  gain 
with  the  significant  hassle  required  to  install,  setup,  and  maintain  this  system. 

2.3.2.  Phase  Diversity  Wave-Front  Sensor 

Possibly  the  simplest  structure  of  all  wave-front  sensors  appears  in  the  phase 
diversity  wave-front  sensor.  This  type  of  sensor  concentrates  on  superior  algorithms  as  it 
is  not  capable  of  the  basic  autocorrelation  algorithms  generally  used  in  wave-front 
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reconstruction  using  the  other  wave-front  sensors  [3].  The  required  maximum-likelihood 
techniques  require  tremendous  processing  power,  as  addressed  on  a  smaller  scale  for  the 
theoretical  sensor,  and  typically  apply  to  offline  de-convolution  of  an  image  rather  than 
real-time  correction  of  wave-front  aberrations  [6]. 

A  beam  splitter  and  a  second  imaging  device  at  a  greater  focal  length  is  all  the 
additional  hardware  required  for  this  sensor  to  estimate  at  least  the  first  2 1  Zemike 
polynomials  [3].  Once  estimated,  the  coefficients  of  the  Zernike  polynomials  allow  for 
de-convolution  of  the  image  with  the  atmosphere,  allowing  for  imaging  when  guide  stars 
are  not  available  [10],  This  process  takes  an  inordinate  amount  of  time,  and  is  not 
capable  of  sustaining  an  adaptive  optics  system  in  real-time  for  fast-moving  objects  or 
rapidly  changing  turbulence;  however,  enough  information  is  available  for  post¬ 
processing  methods.  Low  contrast  scenes  are  still  difficult  to  image  with  this  method  as 
the  SNR  decreases  significantly. 


telescope 
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2.4.  Theoretical  Maximum  Likelihood  WFS 


A  vector-projection  maximum  likelihood  wave-front  sensor  builds  upon  the 
design  of  the  Shack-Hartmann  and  extends  the  SWAT  wave-front  sensor  requiring  no 
major  hardware  changes  from  the  SWAT  design.  This  hardware  setup  provides  the  same 
readout  noise  reduction  as  the  SWAT  sensor,  while  the  algorithm  used  to  detect  tip  and 
tilt  surpasses  centroiding  in  photon  noise  rejection,  particularly  for  extended  objects,  at  a 
cost  of  higher  computation  time  [5]. 

The  hardware  portion  of  this  sensor  adds  an  additional  beam  splitter  just  before 
the  original  Shack-Hartmann  sensor  exactly  as  the  SWAT  wave-front  sensor  does,  with 
the  split  beam  feeding  an  identical,  but  rotated  90°,  array  of  sub-apertures  and  CCD 
elements.  The  CCD  structures  mirror  the  SWAT  device  as  well  by  using  vector  readouts 
of  the  images  creating  projections  of  the  original  image  in  two-dimensions.  The 
algorithm  then  uses  these  projections  independently  for  the  autocorrelation  related 
maximum  likelihood  estimation  of  tip  and  tilt  [5].  Characterization  for  the  setup  and 
some  statistics  already  exist  from  a  previous  work;  therefore,  extension  into  the  tracking 
and  characterization  for  wave-front  sensing  should  be  simpler  [5]. 

2.5.  Comparison  and  Summary 

Limitations  in  current  adaptive  optics  technologies  constrain  the  ability  to 
perform  ad-hoc  imaging  of  fast-moving,  extended,  or  low-contrast  objects.  These 
limitations  generally  stem  from  the  wave-front  sensor,  as  it  is  the  key  component  in  an 
adaptive  optics  system.  Table  1  shows,  using  a  scale  of  Excellent-Good-Marginal-Poor, 
Shack-Hartmann  and  curvature  wave-front  sensors  have  good  perfonnance  for  various 
object  types  and  a  respectable  response  time,  easily  allowing  for  common  use  today  [1,5, 
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16,  17].  The  lateral  shearing  interferometer  and  phase  diversity  wave-front  sensors  have 


other  advantages  as  seen  in  Table  2  that  outweigh  the  detriments  of  such  complex 
systems  [3,  16,  23].  The  theoretical  maximum  likelihood  sensor  provides  excellent 
image  tracking  capabilities  while  maintaining  a  low  complexity  and  high  response  time 
making  it  an  ideal  candidate  for  further  research  [5]. 


Table  1.  Performance  Comparison  for  Common  Wave-Front  Sensors 
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Table  2.  Additional  Known  Advantages  and  Disadvantages  of  Wave-Front  Sensors 


WFS 

Other  Advantages 

Other  Disadvantages 

Shack-Hartmann 

Requires  Small,  High  Contrast 
Object  for  Good  Estimation 

Shearing 

Interferometer 

Very  Adaptable  to  Current 
Environment 

Requires  Extensive  Tuning 

Curvature 

Requires  Tip-Tilt  Estimation 

First  for  Edges  of  Wave-Front 

Phase  Diversity 

Allows  De-Convolution  of  Image 

Only  Offline  Operation 

Maximum 

Likelihood 

Possible  Off-Edge  Lock  Capability; 
Multiple  SW  Realizations  Possible 

Requires  Estimate  of  the  True 
Image 
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2.6.  Possible  Areas  of  Investigation 

Before  research  beings  to  attempt  a  performance  improvement,  a  benchmark  for 
comparison  is  always  a  good  idea.  To  this  end,  a  Cramer-Rao  lower  bound  for  wave- 
front  sensing  should  establish  a  solid  baseline.  Additional  applications  for  the  maximum 
likelihood  wave-front  sensor  are  of  interest,  to  include  integration  with  phase  diversity 
algorithms,  near  and  off-edge  performance  of  guide  stars,  and  multi-spectral  maximum 
likelihood  analysis.  To  make  a  feasible  sensor,  the  algorithm  must  be  capable  of  real¬ 
time  operations  within  a  closed-loop  system  requiring  algorithmic  analysis  and 
decomposition.  Taking  this  decomposition  of  implementable  algorithms  it  should  be 
possible  to  perform  hardware  simulations  and  analysis. 
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III.  Modeling 


Investigating  new  topics  requires  thorough  modeling  of  the  known  environment 
to  guide  the  research  and  provide  adequate  testing  of  results  from  these  analyses  and 
simulations.  This  chapter  of  the  thesis  defines  the  programming  and  simulation  software, 
the  methods  to  generate  realistic  data  within  these  programming  environments,  and  the 
relevant  facts  surrounding  these  modeling  techniques.  Verification  and  validation  for 
expected  perfonnance  of  investigation  results  requires  not  only  the  modeling  capability 
and  understanding  but  also  development  platforms  for  software  and  hardware  simulations 
and  fabrication. 

The  majority  of  software  validation  and  simulation  uses  MATLAB  version 
7.0.4.365  (R14)  Service  Pack  2  with  the  Signal  Processing  Toolbox,  executing  both  the 
simulated  environment  and  sensor  model  under  test.  However,  some  algebraic, 
differential,  and  statistical  validation  uses  Mathematica  version  5.2  for  symbolic 
manipulation  and  verification  of  complex  formulas.  Simulated  hardware  verification 
requires  a  different  development  environment  and  uses  Altera ’s  Quartus  II  version  5.1 
Build  176  for  both  hardware  modeling  and  testbench  simulation.  These  development 
platforms  provide  a  broad  yet  firm  foundation  for  design  and  assessment  of  image  and 
signal  processing  technologies  through  both  software  and  hardware  elaboration  and 
simulation  capabilities. 
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The  objective  in  modeling  images  is  to  provide  the  most  realistic  and  best-case 


scenarios  for  sensor  characterization,  while  providing  the  sensors  with  enough 
information  to  exceed  modern  perfonnance  expectations. 

3.1.  Image  Modeling  Parameters 

As  with  any  modeling,  the  parameters  for  modeling  are  often  more  important  than 
the  modeling  itself,  and  the  parameters  controlling  image  creation  listed  in  Table  3  are  no 
exception. 

3.1.1.  Wavelength 

Since  atmospherically  induced  optical  tilt  bends  different  wavelengths  of  light 
much  like  a  prism,  ideally  a  sensor  should  receive  only  one  wavelength  to  perform 
estimation  as  an  image  further  distorts  when  combining  different  wavelengths.  To  avoid 
further  distortion,  all  created  images  include  the  assumption  that  the  wavelength  is  quasi- 
monochromatic,  including  a  range  of  0.05  pm  of  wavelengths,  and  fixed  both  spatially 
and  temporally. 


Table  3.  Image  Modeling  Parameters  and  Simulated  Ranges 


Parameter 

Description 

Simulated  Range 

Wavelength  of  Light  Received 

Quasi-Monochromatic  and  Fixed 

Sampling 

Nyquist  or  Higher  Sampling  Rate 

1  to  2  times  Nyquist 

Image  Size 

Size  in  Pixels  of  Captured  Image 

8  to  64  Pixels  Square 

Light  Level 

Sum  Total  of  Light  at  Receiver 

Guide  Star:  100  to  1,000  Photons 
Extended  Object:  6,000-20,000  Photons 

Background 

Intensity 

Additive  Stray  Light  in  Receiver 

0  to  1  Photon  per  Pixel 
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While  the  useful  information  contained  in  other  wavelengths  should  produce  similar 
characteristic  results,  the  tilt  information  from  one  additional  wavelength  will  further 
correct  wave-front  error  by  characterizing  the  true  path  of  light  through  the  atmosphere, 
an  effect  not  modeled  or  investigated.  Only  genuine  images  with  real  data  will  define  the 
actual  frequency  of  light  used  for  modeling  as  sampling  requirements  for  real  images 
require  this  infonnation. 

3.1.2.  Sampling 

Once  the  light  passes  through  the  atmosphere  and  enters  the  telescope,  it  is 
necessary  to  sample  the  point  spread  function  (PSF)  appropriately  according  to  the 
Nyquist  sampling  theorem  to  avoid  aliasing  of  frequency  content  in  the  image  [16]. 
Starting  from  the  cutoff  frequency  of  the  optical  transfer  function  (OTF)  of  the  lens 
shown  in  Equation  1,  Nyquist  sampling  chooses  the  minimum  sampling  frequency  to  be 
at  least  twice  this  cutoff  frequency  [16]. 


/«  = 


D 


(1) 


fs  -  2'  fc 


(2) 


where 

fc  =  Telescope  Optic  Cutoff  Frequency  (radians'1) 
fs  =  Sampling  Frequency  (radians'1) 

D  =  Lens  Diameter  (meters) 

X  =  Light  Wavelength  (meters) 
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Given  that  the  wavelength  of  light  remains  constant  for  this  modeling,  any  adjustment 


required  for  Nyquist  sampling  will  occur  either  by  adjusting  the  sampling  rate  or  the  lens 
diameter.  For  an  image  at  a  fixed  distance  from  the  telescope,  Equation  2  combines  with 
Equation  3,  which  assumes  the  small  angle  approximation,  and  then  controls  the 
wavelength  and  aperture  diameter  based  on  the  actual  angular  coverage  of  a  pixel. 


dx  1 
a  =  —  =  — 

Z  fs 


(3) 


where 


a  =  Angular  Coverage  of  Pixel  (radians) 
dx  =  Size  of  Pixel  on  Object  (meters) 
z  =  Distance  to  Object  (meters) 

Over-sampling  has  the  added  benefit  of  aiding  an  interpolator  for  better  estimation 
results,  but  it  also  decreases  the  light  available  to  each  pixel  causing  detrimental  effects 
explained  in  Section  3.5.  Nyquist  sampling  theory  does  not  address  the  resolution  limits 
between  objects  in  the  image;  therefore,  the  Rayleigh,  Dawes,  or  Sparrow  criteria  do  not 
contribute  to  modeling  and  completely  ignored  to  provide  an  idealized  characterization  of 
the  sensors  [15].  Since  the  sampling  frequency  is  twice  the  cutoff  frequency  as 
determined  by  the  diffraction  limited  effect  of  a  telescope  opening,  a  shift  of  one  Nyquist 
pixel  is  analogous  to  a  slope  in  the  frequency  domain  of  n  radians,  or  one-half  of  one 
wave  of  tilt,  via  the  Fourier  shift  theorem  [11]. 

3.1.3.  Image  Size 

Determining  the  actual  image  size  requires  knowledge  of  the  search  area  for 
wave-front  tilt  as  well  as  the  particular  requirements  of  a  wave-front  sensor,  and  in  an 
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effort  to  guide  the  search  area,  a  quick  derivation  of  the  statistical  nature  of  tilt  follows. 


As  often  used  in  a  Monte  Carlo  simulation  when  generating  a  phase  screen,  or  layer  of 
turbulence  causing  tilt  at  a  set  altitude,  the  Cholesky  factorization  of  the  Zemike 
polynomials’  covariance  matrix  multiplies  a  vector  of  zero-mean,  unit  variance  Gaussian 
random  variables  to  create  a  set  of  statistically  accurate  Zernike  coefficients  [16].  The 
tilt  therefore  remains  Gaussian  and  originates  from  the  low  order  elements  of  the 
covariance  matrix,  captured  in  Equation  4  for  Zernikes  two  and  three  [9]. 


<  =  0.4481 


f  D ^  3 


,  {radians 1 


(4) 


where 

a"tih  =  Variance  of  Tilt  (radians') 

D  =  Diameter  of  Aperture  (meters) 
r0  =  Fried  Parameter  (meters) 

This  fonnula  yields  variances,  and  subsequently  standard  deviations,  less  than  one  when 
compared  to  a  wave  of  tilt,  which  is  2n  radians,  for  either  a  large  telescope  or  a  small 
turbulence  coherence  radius  creating  large  values  of  D/r0,  allowing  for  computation  of  a 
search  window.  An  image  supporting  searches  of  plus  or  minus  four  waves  of  tilt  would 
provide  a  worst  case  of  no  less  than  99.99  percent  possible  tilt  coverage,  requiring  a 
search  space  of  plus  or  minus  eight  pixels  [9].  Prior  temporal  analysis  explains  this 
derivation  in  further  detail  and  indicates  that  in  a  closed  loop  system,  as  an  adaptive 
optics  system  provides,  values  of  tilt  beyond  one  to  two  standard  deviations  are 
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exceedingly  unlikely  unless  the  adaptive  optics  system  loses  lock  requiring  a  greater 
search  space  [4]. 

Different  sensors  require  various  image  sizes  to  allow  for  optimal  performance, 
and  although  telescopes  often  include  four-by-four  pixel  images  for  Shack-Hartmann 
sensors  this  image  size  severely  limits  the  range  of  detection  for  tilt  measurements,  thus 
the  current  most  complex  version  of  this  sensor  defines  the  lower  bound  of  an  eight-by¬ 
eight  pixel  image  as  larger  image  sizes  degrade  read-out  performance.  The  theoretical 
Maximum  Likelihood  sensor  relies  on  a  minimum  of  twice  the  number  of  pixels  to 
correlate  with  compared  to  the  desired  search  space  for  extended  objects,  and  to  achieve 
this  sixteen  pixel  search  space,  the  theoretical  sensor  requires  a  minimum  size  of  thirty- 
two  pixels.  With  the  middle  ground  of  image  sizes  fixed,  the  upper  bound  stems  from  a 
forward-looking  perspective  with  respect  to  greater  turbulence  and  superior  accuracy. 

An  important  assumption  as  image  sizes  grow  beyond  approximately  thirty-two  by  thirty- 
two  pixels  is  isoplanism  of  the  observed  wave-front,  which  may  be  true  for  a  natural 
guide  star  (NGS),  is  probably  not  accurate  for  a  satellite  in  orbit,  and  most  likely 
incorrect  for  a  laser  guide  star  (LGS).  It  is  also  important  to  highlight  that  the 
parameterization  of  image  size,  although  somewhat  dependent  upon  sampling,  is 
independent  of  light-level  and  background  intensities. 

3.1.4.  Light-Level  (Total  Intensity) 

Another  parameterized  variable  in  image  creation  is  the  total  intensity  of  the 
image,  or  light-level,  which  the  telescope  controls  based  on  the  object  imaged,  the 
amount  of  light  split  through  the  beam  splitter  to  the  wave-front  sensor,  and  the 
integration  time  of  the  sensor.  To  avoid  temporal  distortions  caused  by  quickly  changing 
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turbulence,  a  short  integration  time  is  desirable;  and  for  modeling  purposes,  all  images 
assume  an  ideal  integration  time  of  100  ps  [16].  Adaptive  optics  entails  obtaining  light 
levels  ranging  from  ten  to  sixty  percent  of  the  total  light  received  by  the  telescope,  with 
the  least  amount  of  light  required  being  the  most  desirable  as  the  light  for  wave-front 
sensing  detracts  from  that  available  to  the  primary  sensor.  To  show  the  performance  of 
all  sensors  in  acceptable  light  levels  and  to  demonstrate  significant  trends  as  light-levels 
increase  or  decrease,  this  modeling  uses  a  modest  search  range  near  the  lowest  light- 
levels  commonly  used.  It  is  interesting  to  note  that,  as  clarified  in  Section  3.5,  when  light 
level  decreases  the  effective  signal-to-noise  ratio  also  decreases  making  a  correct 
estimate  of  the  tilt  less  likely.  This  is  only  one  way  that  the  contrast  ratio,  or  ratio 
between  the  highest  and  lowest  intensities  in  the  image,  changes,  modifying  the 
background  intensity  also  changes  this  hidden  parameter. 

3.1.5.  Background  Intensity 

The  last  considered  parameter  indicates  the  efficiency  of  the  optical  and  electrical 
components  of  an  imaging  system,  and  can  significantly  affect  the  results  of  certain 
sensor  models,  which  expect  a  black  background  to  perform  estimation.  The  background 
light  level  typically  cumulates  from  stray  light  in  the  imaging  system  as  well  as  stray 
electrons  in  the  image  capture  device,  causing  a  lower  contrast  ratio  and  subsequently  a 
lower  SNR.  A  perfect  imaging  system  could  have  a  background  light  intensity  of  zero, 
while  a  very  poor  system  might  aggregate  an  overwhelming  background  of  one  photon 
per  pixel  or  more  for  the  light  levels  in  the  parameterization  range.  This  background 
effectively  resides  beneath  the  signal  represented  by  an  image,  and  can  cause  a  smooth 
function  such  as  a  Gaussian  to  change  shape  significantly. 
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3.2.  Image  Creation 

To  avoid  unknown  effects  in  simulation,  image  creation  requires  explicit 
knowledge  of  the  characteristics  of  an  image  with  nearly  precise  knowledge  of  the 
statistics  expected  from  turbulence  effects  on  that  image  allowing  for  modeling  validation 
and  reliable  simulations. 


3.2.1.  Two-Dimensional  Gaussian  (Simulated  Laser  Guide  Star) 

The  two-dimensional  Gaussian  represented  by  Equation  5  is  possibly  the  simplest 
image  type  to  model,  and  genuinely  represents  the  nature  of  an  artificially  generated  laser 
guide  star,  which  can  allow  wave-front  correction  for  extended  objects.  It  is  important  to 
note  the  assumption  of  independence  and  equality  for  the  variance  in  both  dimensions  of 
the  Gaussian,  which  may  not  be  correct  for  true  hardware  and  atmospheric  turbulence. 
Additionally,  there  is  an  extra  parameter  C  to  adjust  the  light  level  of  the  image,  which 
simply  scales  the  complete  picture  and  does  not  modify  the  Gaussian  in  any  other  respect. 
Figure  10  and  Figure  1 1  illustrate  the  true  fonn  of  this  image  without  noise  for  a  three¬ 


dimensional  view,  projection  images  in  both  axes,  and  a  two-dimensional  representation. 


i(x,y)  =  )  e  2a' 


(5) 


where 


i  =  Representation  of  2-D  Gaussian  Image  (photons) 
x,  y  =  Pixel  Locations  in  the  Image  (pixels,  e  Integers) 

C  =  Total  Intensity  of  Image  (photons) 
a  =  Standard  Deviation  (pixels,  e  Positive  Reals) 

Assuming  a  =  ox  =  ay  and  the  Two  Dimensions  are  Independent 
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Image  of  2-D  Gaussian  (LGS) 


2-D  Gaussian  (LGS) 


Position  {pixels}  Position  {pixels} 


Figure  10.  2-D  Gaussian  and  Vector 
Projections  in  x  and  y  Planes 
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Figure  1 1 .  Observed  Image  of  2-D 
Gaussian  without  Noise 


Aside  from  representing  a  smooth  function  from  which  over-sampling  and  interpolation 
are  simple,  a  further  benefit  of  using  this  model  is  the  ease  of  creating  projections  of  the 
image  discussed  later  in  Sub-Section  3.2.3.  The  Gaussian  has  an  added  advantage  in  that 
it  is  also  a  close  representation  of  a  diffraction  limited  natural  guide  star  once  the  image 
shape  adjusts  to  reflect  the  effects  of  noise,  in  which  case  a  standard  deviation  of  two 
accurately  represents  both  the  NGS  and  LGS  for  modeling  purposes  [4]. 

3.2.2.  Using  Real  Images  or  Real  Data 

Use  of  real  images  requires  more  constraints  than  merely  careful  handling  of 
imprecision  in  the  Fourier  transform,  these  types  of  images  require  proper  centering  to 
provide  fair  statistics,  band-limiting  to  avoid  aliasing,  and  down-sampling  /  up-sampling 
to  meet  sampling  requirements.  Centering  typically  requires  use  of  an  existing  algorithm, 
such  as  the  centroid,  to  adapt  the  image  with  either  Fourier,  or  sine,  interpolation  or 
another  robust  method  such  as  bi-cubic  or  cubic-spline  interpolation.  Although  sine 
interpolation  is  ideal,  this  simulation  uses  MATLAB’s  cubic-spline  interpolation  for  this 
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step  in  modeling  to  prevent  additional  unnecessary  infonnation  appearing  in  the  black 


background  needed  for  down-sampling  and  band-limiting  operations  of  the  image. 

An  important  parameter  that  drives  modeling  is  the  pixel  size  relative  to  the  actual 
size  of  the  simulated  object  as  defined  by  Equations  2  and  3.  This  modeling  investigation 
seeks  to  use  the  Hubble  satellite  for  large  apertures  on  the  order  of  one  meter  and 
subsequently  a  pixel  size  representing  20  cm  [20].  Also  of  interest  is  a  similar  wave- 
front  sensing  problem  in  which  the  aperture  reduces  to  10  cm,  forcing  a  corresponding 
change  in  pixel  size  to  two  meters.  This  pixel  size  determines  the  appropriate  light-level 
for  tracking  or  wave-front  sensing  for  a  real  object  as  summarized  in  Equation  6  [4]. 
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where 


2 

ApiXei  =  Area  that  a  Pixel  Represents  on  Object  (nr/pixel) 

D  =  Diameter  of  the  Aperture  Opening  (m) 
n  =  Number  of  Pixels  in  Array  (pixels) 

Psun  =  Power  of  Sun  at  Earth’s  Surface  ~  1000  (W/m2)  [4] 

A2  =  Bandwidth  of  Light  (Sensor  ~  0.05x1  O'6,  Visible  ~  0.5x1  O'6)  (m)  [5] 
z  =  Distance  to  the  Object  ~  600x10’  (m)  [20] 

At  =  Integration  time  of  Imaging  Device  ~  100x1  O'6  (s)  [5] 
h  =  Planck’s  Constant  ~  6.626xl0"34  (J  s)  [18] 
v  =  Frequency  of  Light  ~  6x10 14  (Hz)  [18] 

R  =  Reflectance  ~  5  to  20  (%)  [18] 

B  =  Light  allocated  by  Beam-Splitter  ~  10  (for  WFS)  (%)  [5] 
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This  equation  involves  several  variables  including  the  ratio  of  the  aperture  and  distance 
to  the  object,  integration  time,  photon  energy,  reflectance,  and  amount  of  light  sent  to  the 
sensor,  most  of  which  are  constant  [4].  This  formula  sets  the  light-level  values  seen  in 
Table  3  for  the  extended  object  scenario,  and  illustrates  the  light  levels  observed  when 
viewing  the  Hubble  space  satellite.  However,  even  with  correct  light-level  calculations,  a 
poor  orientation  of  the  image  may  result  in  low  contrast  for  a  particular  dimension  and  a 
simple  rotation  of  the  image  will  alleviate  a  reduction  in  perfonnance  for  all  models. 

To  avoid  aliasing  while  down-sampling,  convolution  in  space  with  a  sine 
function,  or  multiplication  by  a  two-dimensional  rect  function  in  frequency  representing 
an  ideal  low-pass  fdter,  removes  high  frequencies  beyond  the  down-sampled  image’s 
bandwidth.  If  fdtering  did  not  occur,  higher  frequencies  would  alias  to  lower 
frequencies,  corrupting  the  image  in  the  spatial  domain;  this  aliasing  is  a  type  of  image 
corruption  that  up-sampling  does  not  suffer.  Down-sampling  is  straightforward  for 
Nyquist  sampled  cases;  however,  other  images,  other  samplings,  or  up-sampling  requires 
sub-pixel  information  provided  by  an  interpolator,  and  for  speed  in  modeling  this  step 
uses  MATLAB’s  bi-cubic  interpolation  after  filtering. 

The  optical  transfer  function  of  a  telescope  further  band-limits  the  image 
representing  light  passage  through  the  particular  optic  in  use  and  allowing  for  proper 
diffraction  limiting  effects  caused  by  a  fixed  aperture.  Using  a  standard  diffraction 
limited  OTF,  with  the  factor  to  over-sample  adjusting  the  aperture  diameter  directly,  the 
impulse  response,  or  magnitude-squared  of  the  Fourier  transfonn,  is  the  point  spread 
function,  which  is  a  two-dimensional  Bessel  function,  or  a  perfect  natural  guide  star  [8]. 
Convolution  of  the  image  and  PSF,  or  spatial  frequency  multiplication  of  the  Fourier 
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transform  of  the  image  and  OTF,  provides  a  properly  band-limited  image  with  the 
characteristics  of  the  current  aperture  appropriately  included.  Throughout  the  down- 
sampling  and  band-limiting  processes  to  create  an  image  without  noise  for  simulation,  all 
real  and  complex  data  from  the  Fourier  transforms  passes  through  every  step  to  limit 
errors  due  to  imprecision  in  the  Fourier  transform.  The  images  for  wave-front  sensing 
and  tracking  are  visible  in  Figures  12  through  15. 

Hubble  for  WFS  Hubble  for  WFS 
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Figure  12.  Hubble  for  Wave-Front  Sensing 
&  Projections  in  x  and  y  Planes 


Figure  13.  Observed  Image  of  Hubble  for 
Wave-Front  Sensing  without  Noise 


Hubble  for  Tracking 


Figure  14.  Hubble  for  Tracking  and 
Projections  in  x  and  y  Planes 


Figure  15.  Observed  Image  of  Hubble  for 
Tracking  without  Noise 
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It  is  important  to  note  that  unless  the  ideal  low-pass  filter  is  the  same  or  larger  size  of  an 
up-sampled  PSF,  or  zero-padded  OTF,  by  the  amount  of  down-sampling  required,  there 
are  minor  errors  near  the  edge  of  the  original  down-sampled  image,  which  is  acceptable 
as  long  as  the  errors  are  outside  of  the  window  of  interest  used  for  simulation. 

3.2.3.  Image  Projection  /  Vectorization 

A  projection  of  an  image  is  purely  the  summation  of  an  image  in  one  dimension 
creating  a  vector  representation  of  the  image;  and  furthermore,  this  projection  occurs 
after  any  cropping  of  the  original  image  to  maintain  appropriate  light-levels  [5],  From  a 
hardware  viewpoint,  this  greatly  increases  the  speed  of  image  readout,  which  is  the  main 
limiting  factor  in  the  speed  of  closed-loop  operation;  unfortunately,  this  decreases  the 
light-level  by  one -half  as  discussed  in  Sub-Section  3.5.2  [5]. 

Unrelated  to  noise  statistics,  an  image  projection  implies  independence  between 
the  two  dimensions  of  an  image,  which  is  true  for  operations  limited  to  projections  of  the 
entire  image  in  a  constant  background  as  seen  below  in  Theorem  1.  To  extend  the 
applicability  of  this  theorem,  not  only  images  wholly  contained  in  a  constant  background 
but  also  images  in  a  relatively  low  background  with  minor  fluctuations  exhibit 
dimensional  independence;  however,  extended  objects  do  not  have  dimensional 
independence  as  new  information  enters  and  exits  the  scene.  Although  this  implies  a 
requirement  for  joint  estimation,  the  modeling  here  assumes  dimensional  independence 
as  this  is  theoretically  sufficient  for  simulation  in  a  closed-loop  environment  [5]. 
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Theorem  1 


A  projected  image  is  independent  with  respect  to  the  shift  parameter  in  the 
orthogonal  dimension  as  indicated  by  this  formula: 

oo 

~PX)  =  1 -  Px ,  y  -  Py  )dy 

-co 


where 


i  =  Representation  of  2-D  Image  (photons) 
x,  y  =  Pixel  Locations  in  the  Image  (pixels,  e  Integers) 
fix,  py  =  Shift  in  x  or  y  direction  (pixels,  e  Reals) 

Proof 

This  first  equation  defines  the  starting  point  by  demonstrating  the  two-dimensional 
inverse  Fourier  transform  of  an  image  with  shifts  in  the  x  and  y  directions  [11]. 


a)=  j  ]/(/„/> 


j2PPxfx  +Pyfy  ) f yy) 


I  =  2-D  Fourier  Transfonn  of  Image 
fx,fy  =  Locations  in  Frequency  Domain  (frequency,  e  Integers) 

Next  is  to  integrate  in  the  y  dimension  to  produce  a  projection  in  the  x  dimension. 
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Separation  of  variables  yields  a  smaller  function  integrated  with  respect  to  y. 
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The  integral  on  the  right-hand-side  equals  a  delta  function. 
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CO  CO  CO 

J  -  /w  - 

—co  —co  —co 

Applying  the  sifting  property  of  the  delta  function,  this  selects  fy  =  0  in  the  integral. 

co  co 

J i(x  -  A,  v  -  /?>>  =  | /(/,,0>-'2'AV'2*>C 

-CO  -CO 

Since  the  right-hand- side  is  simply  the  inverse  one-dimensional  Fourier  transfonn  in  the 
fx  direction  with  no  dependence  on  the  shift  in  the  y  direction,  shifts  in  each  dimension 
are  indeed  independent,  which  the  two-dimensional  Gaussian  demonstrates  further  as  its 
one-dimensional  projection  is  simply  a  one-dimensional  Gaussian. 

Q.E.D. 


3.3.  Image  Shifting 

3.3.1.  Shifting  of  a  Known  Function 

The  requirements  for  shifting  a  known  function  restrict  modeling  only  to  the  point 
of  requiring  a  shift  of  the  function  itself,  which  for  a  Gaussian  is  changing  the  mean  of 
the  function.  By  including  the  requirement  of  a  continuous  function,  sub-pixel  shifts, 
which  are  necessary  to  determine  the  true  statistics  of  the  model,  are  also  straightforward. 

3.3.2.  FFT  /  Sine-Interpolation 

For  images  not  generated  from  a  smooth  function,  modeling  requires  another 
method  for  sub-pixel  shifts,  and  the  best  method  for  sub-pixel  shifting  without  knowledge 
of  a  function  is  through  interpolation,  and  for  small  images  sine  interpolation  is  the  ideal 
method  for  accuracy.  A  good  solution  for  implementing  sine  interpolation  is  using  the 
Fourier  transfonn  and  the  relationship  between  shifts  in  the  spatial  domain  and  phase 
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shifts  in  the  spatial  frequency  domain,  as  found  in  almost  any  discrete-time  Fourier 
transfonn  pair  table,  assuming  circular  shifts  are  acceptable  or  use  of  appropriate  zero 
padding  avoids  circular  shifting  [11]. 

3.3.3.  Sub-Pixel  Shift  Step  Size 

Although  simulation  step  size  typically  determines  smoothness  of  results,  this 
modeling  method  requires  different  step  sizes  to  illustrate  different  statistics  representing 
the  performance  of  the  models  appropriately.  As  a  rule  of  thumb,  the  sub-pixel  step  size 
should  be  one-quarter  to  one -tenth  the  pixel  size  as  a  minimum  for  smooth  results  and  no 
smaller  than  the  interpolation  search  size  used  by  the  search  algorithms  and  set  by  the 
CRLB.  Any  smaller  simulation  step  size  would  provide  no  further  insight  beyond 
quantization  error  for  bias  calculations  and  no  insight  beyond  noise  error  as  indicated  by 
the  CRLB  for  noise  calculations. 

3.4.  Calculating  Bias  and  Mean  Absolute  Bias  (MAB) 

The  error  in  the  presence  of  no  noise  is  difficult  to  remove  and  indicates  the  best 
possible  operating  characteristics  of  a  sensor  as  well  as  the  level  of  tuning  required  by  the 
operator  to  achieve  desirable  statistics.  To  calculate  bias,  subtract  the  true  shift  value 
from  the  estimated  shift  value  as  indicated  in  Equation  7;  however,  bias  can  be  deceiving, 
therefore  it  is  better  to  compute  absolute  bias  for  comparison  purposes. 


AbsolueBias  = 


A 


NoNoise 


P 


(V) 


where 

AbsoluteBias  =  The  Absolute  Value  of  the  Error  in  No  Noise  (pixels,  e  Positive  Reals) 
P\ NoNoise  =  The  Estimated  Shift  Value  without  Noise  (pixels,  e  Reals) 
p  =  The  True  Shift  Value  (pixels,  e  Reals) 
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To  reduce  the  complexity  of  results,  averaging  the  x  and  y  dimensions  is  acceptable  as 
long  as  there  is  independence  between  these  dimensions,  otherwise  unexpected  results 
may  surface.  The  bias  should  hold  some  similar  properties  to  the  image,  in  that  if  the 
image  is  symmetric  the  bias  should  be  also,  and  if  the  image  is  higher  contrast,  the  bias 
should  indicate  a  relatively  larger  change  in  some  areas  of  the  curve  if  such  a  change  is 
visible  in  the  search  window. 

The  average  bias  over  a  given  number  of  pixels  is  the  mean  absolute  bias  (MAB) 
and  provides  a  single  number  to  describe  operation  of  the  sensor  for  a  given  region  of  the 
window.  Regions  of  interest  for  the  above  images  include  an  average  over  plus-or-minus 
four  waves  of  tilt  as  this  encompasses  the  entire  window,  and  plus-or-minus  one-half 
wave  of  tilt  as  this  represents  well  over  fifty  percent  of  tilts  seen  in  a  closed  loop  system 
for  reasonable  values  of  D/ro  as  discussed  previously  in  Sub-Section  3.1.3. 

3.5.  Noise  Generation 

3.5.1.  Poisson  and  Bernoulli  Random  Variables 

The  statistics  for  light  are  intuitive  from  the  packet  perspective  of  light  as  each 
photon  interacts  with  objects  such  as  a  beam  splitter  or  charge  couple  device  (CCD)  as  a 
Bernoulli  random  variable  with  a  low  rate  of  success.  One  way  to  approximate  a  Poisson 
random  variable  is  to  sum  many  Bernoulli  trials,  each  with  a  low  rate  of  success,  hence, 
the  overall  statistics  of  light  being  approximately  Poisson  in  nature  [9].  Equation  8  is  the 
general  form  of  a  Poisson  random  variable,  and  is  the  basis  of  the  statistics  required  for 
modeling  of  noise  for  tracking  and  wave-front  sensing.  The  expectation  and  variance 
statistics  for  Poisson  random  variables  indicate  how  the  SNR  increases  as  the  light 
intensity  increases.  Since  the  variance  increases  at  the  same  rate  as  the  mean,  the 
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standard  deviation  increases  as  the  square-root  of  the  mean,  increasing  the  SNR  in  an 


initially  logarithmic,  and  then  nearly  linear  fashion  as  indicated  by  Figure  16.  Over- 
sampling  decreases  the  SNR  on  a  per-pixel  basis  as  the  light  splits  between  more  pixels 
decreasing  the  available  light  per  pixel.  Since  light-level  also  dictates  the  quality  of  a 
captured  image,  the  lowest  SNR  possible  for  a  wave-front  sensor  to  operate  properly  is 
the  ideal  operating  light-level  and  what  this  modeling  attempts  to  parameterize. 

p(d(x, y)\ (*'(*- Px,y-P>)\ Px,Py)j=— — ^ - e-'{*-P*,y-Py)  (8) 


where 


d 

i 


x,y 


P»  fiy 
E[d 1 
VAR[d] 


The  Observed  Intensity  (photons,  e  Integers) 

The  True  Image  Intensities  (photons,  e  Reals) 

Pixel  Locations  in  the  Image  (pixels,  e  Integers) 

Shift  in  x  or  y  direction  (pixels,  e  Reals) 

i  The  True  Image  Intensity  is  the  Mean  (photons,  e  Reals) 

i  The  True  Image  Intensity  is  the  Variance  (photons  ,  e  Positive  Reals) 

SNR  v  Light-Level 
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Figure  16.  SNR  v  Light-Level 
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3.5.2.  Effects  of  Projecting  Images 

With  the  exception  of  the  Shack-Hartmann  sensor  model,  all  sensor  models  use 
vector,  or  projection,  images,  which  cause  some  interesting  effects  for  the  noise  statistics 
to  compare  properly  between  different  sensors.  Building  on  the  assumption  that  each 
pixel  is  independent,  one  can  show  that  the  sum  of  Poisson  random  variables  is  a  new 
Poisson  random  variable,  with  the  new  mean  and  variance  being  the  sum  of  all  means, 
using  either  the  probability  generating  function,  as  shown  in  Theorem  2,  or  indirect 
convolution  and  knowledge  of  Taylor  series  expansions  for  an  exponential. 


Theorem  2 

The  summation  of  Poisson  random  variables,  or  convolution  of  probability  mass 
functions,  is  another  Poisson  random  variable  with  the  new  rate  being  the  sum  of  the 
rates  of  the  summed  random  variables: 


X  d (*’  y)  p(d M  I  (*(*  -  fix )  I P J)  =  p  Z dix’  y)  I  X  ?'(x  _  fix  >  y  -  fir )  I  fix ,  fi 


V  y 


w 


JJ 


Proof 

This  first  equation  is  merely  the  probability  generating  function  redefined  for  the  Poisson 
random  variable  used  in  this  modeling  and  simulation  [9];  note,  Equation  8  defines  all 
parameters  except  for  z  which  is  the  transform  variable. 

Since  a  summation  of  random  variables  is  really  a  convolution,  this  becomes  a  product  in 
the  z-domain  as  indicated  below. 
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y  y 


Now  using  exponential  properties,  the  product  becomes  a  summation  in  the  exponent, 
and  factoring  out  the  (z-1)  term  puts  the  solution  back  into  the  original  form. 


n<w*> = e 


YJi(x-px,y-Py)  (z-l) 


This  summation  indicates  the  new  rate,  mean,  and  variance  are  simply  the  sum  of 
intensities  from  the  original  image,  providing  an  easy  method  of  computing  statistics  for 
projected  image  data. 

Q.E.D. 


In  addition  to  keeping  the  Poisson  statistical  nature  for  a  projected  image  the  light  level 
decreases  by  half,  which  necessarily  decreases  the  SNR  as  defined  in  the  previous  sub¬ 
section.  By  modeling  two  separate  images  at  half  intensity,  including  noise,  then 
producing  projection  images  for  two-dimensions  as  well  as  summing  both  images 
together  to  form  one,  all  sensors  receive  the  same  noise  statistics  while  some  operate  on 
image  projections  and  others  operate  on  a  complete  image.  This  allows  for  accurate 
computation  and  comparison  of  statistics  for  simulation  and  development  purposes. 

3.5.3,  Computing  Noise  Statistics 

There  are  three  main  statistics  for  comparison  between  sensors  when  computing 
with  noisy  data;  however,  they  are  not  unique  and  only  two  of  them  are  useful  to  the 
developer  and  end  user.  The  first  two  statistics  are  nearly  the  same,  as  one  is  simply  the 
square  of  the  other  before  averaging:  mean  absolute  error  (MAE)  and  mean  square  error 
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(MSE).  To  avoid  confusion,  the  MSE  used  for  modeling  is  the  average  square  error  and 
not  the  estimation  technique  familiar  to  researchers  using  signal  processing  estimation 
methods.  To  compute  these  statistics,  refer  to  Equations  9  and  10,  understanding  that  in  a 
similar  manner  to  MAB  these  statistics  are  clearer  when  averaged  over  a  range  of  shift 
values  such  as  one-half  wave  of  tilt  or  four  waves  of  tilt. 


MAE  = 


P\  Noise 


Trials 


N 


(9) 


where 

MAE  =  Mean  Absolute  Error  (pixels,  e  Reals) 

/3\ No^e  =  The  Estimated  Shift  Value  in  Noise  (pixels,  e  Reals) 

ft  =  The  True  Shift  Value  (pixels,  e  Reals) 

N  =  The  Number  of  Trials  (unitless,  Preferred  to  be  a  Power  of  2) 

Ylftyhbe-Pf 

MSE  =  ^ -  (10) 

N 

where 

2 

MSE  =  Mean  Square  Error  (pixels',  e  Positive  Reals) 

Of  the  two  statistics,  MSE  captures  a  broader  view  as  it  is  a  middle  ground  or 
combination  of  the  MAB  and  VAR,  as  indicated  by  Equation  11,  and  is  useful  to  see 
which  of  the  two  statistics  drives  the  resulting  performance  of  the  sensor  [19]. 

MSE  «  VAR  +  (Bias)2  (11) 

where 

2 

VAR  =  Variance  (pixels",  e  Positive  Reals) 

Bias  =  Absolute  Bias  as  Defined  in  Equation  7  (pixels,  e  Positive  Reals) 
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Variance  (VAR)  as  established  in  Equation  12  appears  nearly  the  same  as  MSE 


with  two  significant  differences:  1)  the  sample  mean  is  the  subtrahend  rather  than  the 
true  shift  value,  and  2)  the  divisor  after  summing  the  sample  is  one  less  than  the  total 
number  of  trials  making  it  an  unbiased  estimate  of  the  variance.  As  a  second  order 
statistic,  VAR  indicates  how  well  a  sensor  can  perfonn  for  a  given  noisy  environment, 
and  is  impossible  to  remove  without  changing  the  type  of  estimation  or  optical  setup. 


£ 

Trialsl 


VAR  = 


fi\Noise 


Z4 


Noise 


Trials 


N 


N- 1 


(12) 


Although  only  qualitative  bounds  are  available  for  average  bias  and  error,  it  is  possible  to 
provide  an  analytical  bound  for  variance  that  defines  the  efficiency  of  an  algorithm’s 
ability  to  reject  noise  in  various  conditions.  With  the  proper  background  and  modeling 
capabilities,  this  Cramer-Rao  lower  bound  can  provide  insight  into  development  of  an 
algorithm  to  improve  tracking  and  wave-front  sensing,  while  verifying  simulation  and 
experimental  results. 

3.6.  Summary 

Accurate  modeling  not  only  guides  research  to  feasible  solutions  but  also 
provides  a  method  to  verify  research  results  before  actual  implementation.  For  this 
research  effort,  generation  of  images  and  the  noise  statistics  that  surround  them  is  the  key 
to  better  understanding  and  estimation  of  wave-front  parameters  and  tracking  shifts. 
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IV.  Analysis 


Investigation  in  wave-front  sensing  requires  thorough  knowledge  of  the  current 
estimation  techniques,  environmental  parameters,  and  modeling  practices  to  provide 
guidance,  insight,  and  validation  capabilities  for  research.  The  key  areas  for 
investigation  for  this  research  include  bounds  on  variance  to  quantify  performance  for 
any  wave-front  sensor,  search  algorithm  optimization  for  the  maximum-likelihood  wave- 
front  sensor  to  meet  or  exceed  timing  requirements,  and  an  implementation  proposal  with 
hardware  realization  of  the  sensor  algorithm  to  demonstrate  feasibility  of  this 
implementation.  Theory  can  provide  excellent  guidance  for  algorithm  development  and 
hardware  implementation  if  applied  correctly,  as  this  research  attempts  to  do;  and  the 
proper  use  of  theoretical  results  can  significantly  shorten  development  time  compared  to 
trail  and  error  analysis. 

4.1.  Cramer-Rao  Lower  Bound  (CRLB)  for  Tilt  Estimates  Obtained  with  LGS 

4.1.1.  Relevant  Statistics,  Assumptions,  and  Setup 

As  the  CRLB  is  a  bound  on  variance,  it  requires  statistical  background  and  noise 
information  given  a  particular  type  of  data  and  a  proper  foundation  to  provide  meaningful 
information.  The  basis  of  analysis  resides  with  Equation  8  in  Chapter  III  and  the 
assumption  of  a  form  of  the  laser  guide  star  for  the  image  as  a  two-dimensional  Gaussian 
represented  by  Equation  5;  repeating  both  equations  below  provides  clarity. 
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(8) 


d  =  The  Observed  Intensity  (photons,  e  Integers) 
i  =  The  True  Image  Intensities  (photons,  e  Reals) 
x,  y  =  Pixel  Locations  in  the  Image  (pixels,  e  Integers) 
fix,  fiy  =  Shift  in  x  or  y  direction  (pixels,  e  Reals) 

-(py) 

i(x,  y)  =  c(2;rcr2 )  '  e  2o“ 

where 


(5) 


i  =  Representation  of  2-D  Gaussian  Image  (photons) 
x,  y  =  Pixel  Locations  in  the  Image  (pixels,  e  Integers) 

C  =  Total  Intensity  of  Image  (photons) 
a  =  Standard  Deviation  (pixels,  e  Positive  Reals) 

Assuming  a  =  ax  =  ov  and  the  Two  Dimensions  are  Independent 
Assuming  that  the  dimensions  are  independent,  using  the  fact  that  the  sum  of  Poisson 
random  variables  is  another  Poisson  random  variable,  and  using  the  image  projection 
technique  to  remove  one  of  the  dimensions,  Equations  8  and  5  become  marginal  with 
respect  to  x  in  Equations  13  and  14  below. 

p(d(x) I  (;(*  -  A )  |  A ))  =  ]  ( 1 3 ) 

1  Z^. 

i(x)  =  C^Ikcj]  )  2  e 2tJ'  (1 4) 
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Assuming  the  pixels  for  a  projected  image  are  independent  this  results  in  the  joint 
probability  mass  function  (PMF)  representing  the  joint  a  priori  density  function  in 
Equation  15. 

II M40 1  Oh  -  A )  i  A )) = 1  (i5> 

The  random  variables  in  this  equation  are  the  shift  represented  by  [3X  and  two  unwanted 
parameters  C  and  Oi,  which  are  part  of  the  assumed  image  fonn.  However,  this  equation 
does  not  account  for  windowing  of  either  the  data  or  initial  image  as  the  information 
captured  is  finite  in  size,  and  the  proposed  maximum-likelihood  sensor  requires  further 
windowing  to  search  over  shifts  and  to  prevent  detrimental  effects  from  new  data 
entering  the  scene  [5].  To  limit  the  product  properly,  a  windowing  function  on  both  the 
true  image  as  well  as  the  captured  image  simply  bounds  the  limits  for  the  product 
function,  and  completes  the  probability  information  required  to  derive  the  CRLB  for  an 
unbiased  estimator. 

4.1.2.  Derivation 

As  noted  previously,  estimation  of  the  shift  parameter  is  the  goal;  however,  two 
additional  parameters  require  estimation  as  well  and  therefore  a  joint  estimation  approach 
of  these  parameters  and  the  CRLB  serves  as  an  accurate  lower  bound  for  Gaussian 
images.  To  derive  a  CRLB  requires  computation  of  the  elements  that  compose  the  Fisher 
Information  Matrix  as  defined  in  Equation  16,  where  the  diagonal  elements  of  the  inverse 
of  this  matrix  are  the  CRLB  for  the  respective  parameters  in  Equation  17  [19]. 
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(16) 


j  dApAj 

where 

J  =  Elements  of  the  Fisher  Information  Matrix 
D  =  the  Entire  Vector  d(x) 

A  =  the  Parameters  to  Estimate  (fix-,  <Ji,  C) 

(17) 

where 

J  =  The  Fisher  Information  Matrix 

As  this  equation  calls  for  the  log  of  the  joint  likelihood,  Equation  18  illustrates  the  log  of 
Equation  15,  with  further  simplifications. 

1 1"0V W I  Of  -  A )  I  A )))  = 

=  X  (ln(/(x  -  fix  )rf(jr))-  ln(j(x)!)  +  ln(e_'(*"A))) 

= X  M  ln(M  _  fi  ))  ”  lrM MO  -  M  -  fix ))  ( 1 8) 

To  reflect  the  additional  unwanted  parameters,  Equation  19  includes  a,  and  C  as 
additional  givens  in  the  log-likelihood,  where  the  true  image  conditioned  on  these 
parameters  represent  the  vector  A  in  the  Fisher  Information  Matrix  and  the  observed 
data  represents  D . 

X  M/M M I  (*(*  -  fix )  I  fix, i c))) = X (d M|n(/(*  -  fix ))  -  |nM MO -  i(x  -  fix ))  ( i 9) 
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An  assumption  that  windowing  equation  19  does  not  change  the  derivative  allows 


computation  of  partials  without  knowing  the  derivative  of  the  window  function;  however, 
this  assumption  is  only  a  close  approximation  when  the  observed  image’s  intensity 
decreases  to  zero  at  the  edge  of  the  window  making  the  CRLB  applicable  for  images  with 
a  large  black  background  or  zero-shift  estimation.  This  is  the  best-case  operation  of  a 
sensor,  and  still  provides  an  accurate  lower  bound  for  performance  of  estimation 
techniques.  As  the  partial  derivatives,  logarithm,  and  partial  derivatives  of  the  log  of  the 
Gaussian  image  form  appear  several  times  in  the  next  derivation,  Equations  20  through 
26  summarize  these  results  based  on  Equation  14. 


a  ,•(*-&)=. 6- a  )-(*~A) 


dp 


(20) 


o, 


do 


o, 


O  : 


=  i(X~Px) 


{x~PxY  ~<r? 


(21) 


-jL'ix-Px)  =i(x-Px)y, 

dC  C 


(22) 


ln(/'(x  -  Px  ))  =  Inf  dlTiof) 


i  A 


(x  ~  P,  )2 

2cr,2 


(23) 


4—\n{i{x-pP) 

dpx 

-P—  ln(/(x  -  Px  )) 

do, 

jUn(/(x-/?J) 


(x~Px) 

(x-Px)2  1  _  (x  -  P,  )2  - 

<t3  '  <r3 

i 

c 


(24) 


(25) 

(26) 
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Leveraging  the  information  in  Equations  20  through  26,  the  following  equations  compute 
the  first  partials  of  the  log-likelihood,  which  are  also  useful  for  maximum-likelihood 
estimation  of  these  parameters  with  the  Gaussian  image  assumption. 

5 


-  X  In (p(d (x)  I  (i(x  -  px )  I  px ,  cr . ,  C ))) 


0A 


-  X  (d  (x)ln(z(x  -  fl ))  -  In (d (x))  -  i(x  -  fl )) 


^-(J(x)ln(?'(x  -  A )))“  T^-ln(^(x))-  -^-i(x  -  A ) 

,v  l^A  5A  3A 


cr. 


^  y 


(27) 


5cr, 


X  ln(p(rf  (x)  |  (i(x  -0x)\  px  ,o„C ))) 


do. 


-  X  (rf (x)ln(i(x  -  P, ))  -  In (d (x)!)  -  i(x  -  px )) 


X  [ d(x) ln0(x  -  Px  )) -  -p—  ln(d  WO  -  p—  *(x  -  Px ) 

,  ^dcr,.  do do  I  J 


-lU> 

cr,  cr, 


J 


(28) 


dc 


X In (p(d (x)  I  (i(x  - Px )  I  px ,o,,C ))) 


dc 


X  (d (x)ln(z(x  -  px ))  -  In (d (x)) -  i(x  -  px )) 


Z  f  d Wln(*(x  -  A ))  -  ^  ln(^  WO  ~  A ) 


48 


Unfortunately,  the  Fisher  Infonnation  Matrix  requires  the  second  partials  with  respect  to 
all  estimated  parameters;  therefore,  it  is  a  square  matrix  and  the  elements  of  the  matrix 
should  be  symmetric  about  the  diagonal  as  the  order  of  the  partial  derivatives  should  be 
reversible. 
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Substituting  Ufa  )  =  -  -L 
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Because  the  second  partial  derivatives  are  interchangeable  in  order,  this  will  create  a 
symmetric  Fisher  Information  Matrix  as  expected,  further  corroborating  the  results 
above. 
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The  final  step  to  complete  the  elements  of  the  Fisher  Information  Matrix  is  to  take  the 
negative  expectation  of  each  second  partial  derivative,  recognizing  that  the  only  random 
variable  in  Equation  19  is  d(x),  whose  expectation  is  Equation  39. 
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To  put  the  result  in  Equation  40  in  perspective,  if  the  parameters  a,  and  C  are  given,  then 
the  inverse  of  this  would  be  the  CRLB  for  the  single  parameter  estimation;  however,  joint 
estimation  requires  the  rest  of  the  terms  as  well. 
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There  are  no  random  variables;  therefore,  the  expectation  has  no  effect. 
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There  are  no  random  variables;  therefore,  the  expectation  has  no  effect. 
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To  summarize  these  results,  Equation  46  displays  the  entire  Fisher  Information 
Matrix,  and  as  Equation  17  illustrates,  each  diagonal  element  of  the  inverse  of  this  matrix 
is  the  CRLB  for  the  respective  estimated  parameters.  The  off-diagonal  elements  of  the 
inverse  Fisher  Information  Matrix  represent  the  bounds  on  covariance  terms  determining 
independence,  or  lack  of  independence,  between  the  estimated  parameters. 
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This  matrix  is  too  complex  to  take  the  inverse  of  symbolically;  however,  assuming 
independence  between  parameters,  as  some  off-diagonal  elements’  anti- symmetric  nature 
indicates,  may  allow  easy  inversion  of  the  diagonal  elements.  Figure  17  in  the  next  sub¬ 
section  illustrates  the  numerically  calculated  results  for  fixed  parameters.  Computing 
results  numerically  indicates  two  main  points:  for  aliased  images,  small  images,  and  near 
the  edge  of  a  fixed  window,  the  bound  appears  incorrect;  and  a  point  solution  for  a  zero- 
shift  estimate  appears  valid  for  the  majority  of  the  window. 

4.1.3.  Simplification  and  Further  Assumptions 

It  is  possible  to  compute  a  zero-shift  solution  for  the  CRLB  as  the  model  is 
accurate  for  this  condition  since  a  zero-shift  produces  minimal  discontinuities  due  to 
windowing  in  the  derivative.  Since  the  sampling  of  the  Gaussian  shape  produces  nearly 
linear  regions  between  each  sample,  the  summations  over  the  values  of  x  emulate 
integrals. 
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Select  the  following  variables  for  integration  by  parts  (IBP)  and  pull  the  constants 
out  in  front  of  the  integral. 
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Attempting  to  evaluate  the  integral  for  U  V  yields  oo/oo;  therefore,  L’Hopital’s 
Rule  can  still  provide  the  limit  as  this  function  approaches  oo  in  both  directions. 
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Rearranging  this  result  produces  a  constant  multiplied  by  the  integral  of  a 
Gaussian,  which  is  just  the  constant. 
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Selecting  the  following  variables  for  IBP  for  the  left-most  integral  creates  a 
positive  two  times  the  integral  on  the  right  after  one  step  of  IBP,  which  combine  to 
produce  a  single  positive  integral  as  shown  below. 
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Rather  than  laboriously  apply  L’Hopital’s  Rule  three  times  to  determine  the  limit 
of  this  fraction  as  it  approaches  oo  in  both  directions,  it  is  clear  that  the  numerator  will 
eventually  be  a  constant  and  the  denominator  will  remain  an  exponential  dependent  on  x, 
again  yielding  0-0. 
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The  right-hand  integral  is  simply  an  integral  of  a  constant  multiplied  by  a 
Gaussian,  which  is  the  constant,  whereas  the  center  integral  is  identical  to  Jn  and 
therefore  equal  to  C  o;  .  The  left-hand  integral  requires  temporarily  pulling  the  constants 
out  front  and  integration  by  parts  as  shown  below. 
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Again,  L’Hopital’s  Rule  indicates  the  first  term  approaches  zero;  however,  re¬ 
arranging  the  remaining  integral  reveals  the  same  form  as  Jn,  again  simplifying  the 
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integration  process  by  providing  the  answer  of  C  of. 
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Inserting  these  results  into  the  original  equation  yields  the  following: 
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The  left-most  integral  is  identical  to  the  form  found  in  Jn;  therefore,  using  the 
same  solution  and  L’Hopital’s  Rule  calculation  results  in  the  following: 


Equation  53  summarizes  these  results  for  the  Fisher  Information  Matrix  indicating  that 
the  magnitude  of  the  true  Fisher  Information  Matrix  is  less  than  or  equal  to  this 
approximation  to  provide  a  true  lower  bound  and  that  the  parameters  are  uncorrelated, 
with  the  matrix  inverse  shown  in  Equation  54  providing  an  accurate  approximation  of  the 
CRFB  for  all  three  parameters. 
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Figure  17  below  indicates  the  match  of  these  zero-shift  solutions  to  the  original 
numerically  computed  CRLB  from  Equation  46,  which  verifies  the  theoretical  results, 
while  Appendix  A  contains  a  Mathematica  notebook  to  ensure  every  step  is  correct. 

CRLB  v  Shift 


Figure  17.  CRLB  Numerical  and  Analytical  Solution 
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Additionally,  the  last  estimated  parameter,  C,  directly  stems  from  the  variance  of  a 
Poisson  random  variable  as  summing  all  intensities  in  the  vector  produces  this  as  the 
variance  of  the  new  Poisson  random  variable  as  described  in  the  Section  3.5. 

4.1.4.  Benefits  and  Discussion 

The  above  theoretical  bound  on  variance  for  shift  estimation  has  two  main  areas 
of  benefit  due  to  the  simplicity  and  completeness  of  the  bound;  first,  the  bound  can  guide 
researchers  in  implementation  of  estimation  algorithms,  and  second,  it  can  guide 
technicians  towards  reasonable  light-levels  and  images  for  shift  estimation.  The 
implementation  benefits  are  two-fold  in  that  estimation  algorithms  that  search  for  sub¬ 
pixel  shifts  need  only  search  to  the  square  root  of  the  minimum  variance  given  by  the 
bound  as  noise  error  overrides  any  quantization  error  in  the  model.  The  bound  also 
provides  an  analytical  method  to  validate  the  sensor  model  and  simulation  results  by 
determining  if  the  model  is  efficient  in  achieving  the  bound  and  providing  another  fonn 
of  verification  for  modeling  by  allowing  comparisons  to  this  independent  bound. 

4.2.  Maximum  Likelihood  Optimized  Search  Algorithm 

4.2.1.  Relevant  Statistics,  Assumptions,  and  Setup 

Leveraging  the  noise  statistics  from  Chapter  III  and  the  projection  of  an  image 
exhibiting  these  statistics  at  the  beginning  of  this  chapter,  the  maximum-likelihood  (ML) 
estimator  uses  the  joint  a  priori  distribution  as  shown  in  Equation  15  to  determine  what 
the  estimate  should  be  according  to  the  criterion  in  Equation  55. 

A  ml  =  ar§  max  Ff  P(d  M I  (*(*  ~PX)\  Px ))  (55) 
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To  minimize  the  computational  complexity  of  performing  numerous  multiplications 


during  a  search  over  values  of  px,  the  same  result  is  available  from  the  natural  log  of 
Equation  55  as  indicated  in  Equation  56. 

PxML  =  ar8  max  X  Xn(p(d  W  I  (i(X  -Px)\  Px  )))  (56) 


This  maximum-likelihood  approach  differs  from  the  maximum  a  posteriori  (MAP) 
estimation  approach,  which  seeks  to  maximize  the  likelihood  over  the  joint  a  posteriori 
distribution  found  in  Equation  57  by  using  either  Bayes’  Rule  or  applying  the  Law  of 
Total  Probability  and  the  criterion  found  in  Equation  58  [5]. 


n  />((*(*  -  p*)\  Px)  i  d(x\p'x)= n 


p(d (x) I  (/(x - px ) I  Px ))p((i(x -px)\px)\0x) 
p(d{x)) 


(57) 


Pxmap  =  argma xj] />((*(*  " Px)  I  Px)  I  d(x\Px)  (58) 

x 


Since  this  technique  seeks  to  maximize  the  joint  a  posteriori  distribution  for  a  given 
value  of  px,  the  marginal  with  respect  to  the  observed  image  is  unnecessary,  while  using 
the  log  of  the  a  posteriori  distribution  and  expanding  further  simplifies  the  search  as 
indicated  by  Equation  59. 

Pxmap  =  ar8  max  X ln^ M I  (?'(*  "  Px  )  I  Px  )))  +  X  -  Px  )  I  Px  )  I  P'x  ))  (59) 


In  the  case  that  the  right-hand  term,  or  the  prior  probability  of  px  given  the  previous  shift, 
is  uniform,  the  dependence  on  the  prior  withdraws  causing  the  MAP  and  ML  estimates  to 
become  equal.  As  mentioned  in  the  chapter  on  modeling,  this  prior  distribution  is 
Gaussian  in  nature;  however,  the  parameters  required  for  this  distribution  are  not 
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available  leaving  the  choice  of  assuming  a  uniform  distribution,  as  previously  performed 
in  literature  [4,  5]. 

The  expansion  of  Equation  56  yields  the  log-likelihood  for  implementation  and 
the  first  part  of  the  MAP  estimator,  should  the  prior  information  become  available,  as 
specified  by  Equation  60. 


P 


XML 


arg  max 


1 1« 


d{x)\ 


-ifi-A) 


=  arg  max  ^  (ln(/(x  -  fix  )'/,',)+  ln(e  >ix  A  * ) -  1  n (<;/ (x))) 

X 


=  arg  max  ^  (d (x)ln(/(x  -  J3X ))  -  i(x  -  ft x )  -  In (7/ (x))) 

X 


Since  the  final  term  does  not  depend  on  px,  the  tenn  drops  out. 

=  arg  max  £  (d (x) ln(/(x  -  J3X ))  -  i(x  -  P x ))  (60) 

a: 


Ideally,  the  most  efficient  method  for  obtaining  the  shift  estimate  is  through 
taking  the  derivative  of  the  above  function,  setting  it  to  zero,  solving  for  the  nodes  of  the 
function,  and  detennining  which  node  has  the  largest  peak.  It  may  be  possible  to  derive  a 
closed-form  solution  for  the  derivative  of  the  log-likelihood  provided  the  derivative  of 
the  original  image  i(x-[3x)  also  has  a  closed  fonn  solution,  which  is  both  image  and  shift 
specific.  This  typically  is  not  possible;  however,  leveraging  the  unique  properties  of  the 
Fourier  Transform  and  the  interchangeability  of  derivatives  and  integrals  with  additional 
assumptions  may  provide  such  a  closed  form  solution  for  faster  estimation. 

There  are  numerous  alternative  approaches  to  searching  the  log-likelihood 
including  using  a  known  function  for  the  true  image  to  allow  use  of  all  of  the  data  in  the 
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log-likelihood  calculation.  Although  a  known  function  could  preclude  the  need  to 
window  the  data,  simply  adjusting  the  algorithm  to  accommodate  both  a  fixed  true  image 
size  and  a  dynamic  true  image  size  allows  a  greater  search  area  for  smaller  image  sizes 
and  possibly  better  performance  without  modification  to  most  software  or  hardware.  An 
additional  assumption  when  working  with  actual  data  is  regarding  which  projected  image 
to  shift  for  searching  over  different  search  estimates,  as  it  is  possible  to  shift  either  the 
true  image  or  the  observed  data.  Theory  indicates  that  the  random  parameter  is  in  the 
original  image,  and  therefore  shifting  of  the  true  image  is  appropriate;  however,  it  may  be 
interesting  to  characterize  the  noise  rejection  capabilities  of  shifting  the  observed  data 
also,  as  hardware  interpolation  is  possible  for  this  sub-pixel  shift  technique.  This 
research  focuses  on  implementing  a  maximum-likelihood  search  approach  using  the  joint 
log-likelihood  defined  by  Equation  60  in  an  efficient  manner  to  temporally  compete  with 
the  Shack-Hartmann  and  SWAT  wave-front  sensors,  as  current  research  indicates  a 
statistical  performance  improvement  with  the  ML  sensor  for  extended  objects  [5]. 

4.2.2.  Properties  of  Log-Likelihood  Leveraged 

As  observed  from  computing  sample  log-likelihoods  using  the  modeling 
techniques  in  Chapter  III,  there  are  several  properties  of  the  log-likelihood  curve  that 
lend  themselves  to  an  optimized  search  algorithm.  Figure  18  and  Figure  19  illustrate  the 
log-likelihood  curve  for  the  laser  guide  star  as  described  in  the  modeling  chapter.  The 
most  important  attributes  include  the  large  main  node,  significantly  smaller  nodes  and 
distortion  due  to  noise,  and  mild  peaks  at  the  end-points  indicating  perfonnance  in  a 
constant  background. 
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Log-Likelihood  for  LGS 


Log-Likelihood  for  LGS 
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Figure  18.  Log-Likelihood  for  Gaussian  in 
Noise 


Figure  19.  Log-Likelihood  for  Shifted 
Gaussian  in  Background  and  Noise 


The  peak  of  the  main  node  is  the  estimate,  and  a  search  algorithm  that  rejects  the  noise 
and  other  characteristics  of  the  log-likelihood  could  search  in  an  efficient  manner  using 
the  concave-down  properties  of  the  main  node.  The  peaks  at  the  edge  of  the  log- 
likelihood  window  develop  when  the  black  background  of  the  true  image  covers  greater 
than  one -half  of  the  search  window  making  it  more  likely  that  the  object  has  moved 
completely  out  of  view.  These  end-point  peaks  have  the  unique  characteristic  of  being 
slightly  greater  on  the  side  of  the  log-likelihood  curve  that  contains  the  main  node  for  a 
significant  shift,  also  lending  possible  simplification  to  the  search  algorithm. 

4.2.3.  Search  Algorithm  Definition 

The  goal  of  this  search  algorithm  is  to  perform  an  efficient  search  of  the  concave- 
down  portion  of  the  main  node  of  the  log-likelihood,  while  being  robust  enough  to  reject 
interference  from  noise  and  other  artifacts  in  the  search  window.  A  method  that  skips  the 
noise  and  artifacts  by  quickly  finding  the  main  node  before  performing  finely  stepped 
search  operations  effectively  meets  these  requirements  and  provides  a  robust  solution  for 
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the  search  algorithm.  By  pre-computing  the  interpolated  true  image  and  its  logarithm, 
each  sensor  can  use  this  data  without  wasting  more  computations  allowing  the  optimized 
algorithm  to  focus  on  searching  the  log-likelihood  and  ignoring  the  requirements  of  the 
input  data.  The  implemented  search  algorithm  has  two  phases;  the  first  grid  search  is 
optional  depending  on  the  shape  of  the  log-likelihood  and  the  second  implements  an 
optimal  search  algorithm  requiring  minimal  memory  storage  for  a  known  concave  down 
function  as  indicated  by  the  program  flow  in  Figure  20. 

There  are  two  different  ways  to  describe  the  main  search  algorithm  using  modern 
search  techniques;  the  first  method  stems  from  the  Gradient  Decent  algorithm,  while  the 
second  method  builds  upon  the  algorithmic  concept  of  a  Binary  Search  Tree.  From  the 
Gradient  Decent  perspective,  which  is  the  basis  for  development,  this  algorithm  perfonns 
Gradient  Ascent  by  climbing  the  log-likelihood  curve,  where  the  step  size  and  slope 
determinations  are  the  unique  and  key  components  of  the  algorithm.  The  step  size  uses 
Bisection  to  determine  the  next  point  in  the  search,  as  it  is  easy  to  compute  this 
dynamically  changing  step  size  and  reduces  the  complexity  of  the  search  significantly. 
The  slope  detennination  ensures  locating  the  peak  by  determining  which  direction  to 
climb  when  encountering  a  larger  log-likelihood  value  with  version  2  of  the  Select 
Window  Endpoints  block  choosing  the  new  search  region.  The  algorithm  properly 
assumes  that  the  slope  is  toward  the  current  largest  value  if  the  new  log-likelihood  of  the 
computed  point  is  less  than  the  current  maximum  and  performs  endpoint  detection  using 
version  3  of  the  Select  Window  Endpoints  block. 
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Figure  20.  Flow  Diagram  of  Optimized  Log-Likelihood  Search  Algorithm 
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This  slope  determination  allows  the  search  to  narrow  the  search  region  for  the  peak 
quickly  and  completes  the  Gradient  Decent  algorithm.  From  the  perspective  of  dividing 
the  problem  into  search  regions,  another  view  of  the  algorithm  emerges  in  the  form  of  a 
Binary  Search  Tree,  with  the  first  node  being  the  entire  search  window,  this  node’s 
children  being  the  left  and  right  halves  of  the  search  window,  and  continuing  until  only 
individual  elements  are  the  leaves  of  the  tree.  As  the  search  progresses  the  algorithm 
makes  a  decision  at  each  node  to  guide  which  children  to  select  and  proceeds  with  a 
depth-first  search  of  the  entire  tree;  and  since  these  decisions  are  final,  upon  reaching  a 
leaf,  the  index  of  the  leaf  is  the  result  of  the  search  algorithm.  This  search  algorithm  has 
the  added  advantage  that  it  requires  memory  for  only  three  index  values,  log-likelihood 
values,  and  their  attributes  making  it  very  feasible  for  implementation  in  a  compact 
embedded  architecture  for  fast  operation. 

The  optional  grid  search  allows  the  algorithm  to  proceed  given  the  main  node  of 
the  log-likelihood  cannot  be  found  within  the  first  three  computations  of  the  bisection- 
based  search  algorithm.  This  search  is  simply  a  linear  search  across  the  window,  with  the 
results  guaranteeing  a  narrowed  search  region  decided  by  version  1  of  the  Select  Window 
Endpoints  block  including  only  the  main  node  of  the  log-likelihood  curve.  The  nature  of 
this  search  requires  one  additional  memory  location  to  ensure  proper  capture  of  the 
maximum  log-likelihood  and  the  window  surrounding  this  maximum,  but  it  remains 
simple  due  to  capturing  the  maximum  while  computing  the  new  log-likelihood  values 
and  performing  only  one  slope  detection  when  complete. 

These  algorithms  require  knowledge  of  the  size  of  the  observed  image,  or  data, 
and  true  image  arrays  as  well  as  the  number  of  sub-pixel  points  interpolated  between  the 
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points  in  the  real  data;  however,  they  do  not  require  knowledge  of  the  type  of  interpolator 
used,  which  Appendix  B  further  confirms.  The  number  of  sub-pixel  points  between  real 
points  is  a  direct  result  from  the  CRLB  analysis,  as  selecting  an  interpolation  allowing 
searches  smaller  than  the  minimum  standard  deviation  simply  increases  computation  time 
with  no  additional  accuracy  benefits.  Using  twice  the  light-levels  suggested  in  the 
chapter  on  modeling  as  well  as  the  assumption  of  a  standard  deviation  of  two  for  a 
Gaussian  image  results  in  a  minimum  standard  deviation  of  error  for  shift  estimation  on 
the  order  of  0.03  pixels.  Interpolating  by  25  points  narrowly  accommodates  this  value, 
while  26  sufficiently  surpasses  this  minimum  error;  therefore,  a  step  size  of  the  inverse  of 
27  represents  an  ideal  interpolation  level  to  allow  for  proper  estimation  of  any  type  of 
image  while  only  requiring  one  additional  log-likelihood  computation.  Since  each  sensor 
requires  only  one  true  image  as  indicated  in  Section  4.3,  a  separate  processor  or  module 
optimized  for  interpolation  could  provide  this  information.  Research  indicates  that  linear 
interpolation  is  good  enough  for  estimation;  however,  a  characterization  of  the 
performance  difference  between  linear  and  the  best  possible  cubic-spline  interpolator  that 
MATLAB  includes  is  an  interesting  area  of  investigation  included  in  this  research  [5]. 

4.2.4.  Computational  Complexity  Analysis  and  Comparison 

Computational  complexity  research  provides  insight  into  how  difficult  or  how 
long  an  algorithm  takes  compared  to  another  algorithm  given  the  same  inputs,  and  the 
following  analysis  performs  this  comparison  for  the  Shack-Hartmann,  SWAT,  non- 
optimized  Maximum-Likelihood,  and  optimized  Maximum-Likelihood  wave-front 
sensors.  This  is  a  comparison  of  the  software  portion  of  the  algorithm  only;  the  delays 
associated  with  hardware  CCD  readout  and  data  transfer  do  not  appear  in  these 
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computations  or  comparisons  and  are  important  considerations  outside  of  this  discussion. 


The  standard  notation  of  On(n)  describes  the  computational  complexity  by  indicating  the 
asymptotic  nature  of  a  function  based  on  the  subscripted  variable.  To  clarify  the 
calculations,  each  sensor’s  computational  complexity  stems  from  a  single-dimension 
estimate  of  the  wave-front  tilt,  with  significant  constants  indicated  in  the  standard 
notation  with  a  for  multiplication  or  ‘+’  for  addition  between  the  two  independent 
complexities.  The  Shack-Hartmann  wave-front  sensor  requires  software  projection  of  the 
image  when  computing  the  centroid;  therefore,  this  excess  computation  time  appears  in 
the  number  of  additions  required  for  an  image  as  indicated  by  Table  4.  The  faster  SWAT 
wave-front  sensor  capitalizes  on  hardware  projection  of  the  image,  which  also  decreases 
the  read-out  time  of  the  data,  and  shows  a  linear  growth  of  both  additions  and 
multiplications  with  image  size.  The  far  more  complex  maximum  likelihood  algorithm 
breaks  into  three  phases,  with  one  optional  phase,  where  the  summation  of  these  three 
phases  indicates  the  total  complexity  of  the  search  algorithm.  Assuming  a  fixed  window 
of  one-half  the  image  length,  the  complexity  of  the  pre-compute  phase  and  a  non- 
optimized  search  algorithm  appears  in  Table  5,  while  Table  6  presents  the  complexity  of 
the  optimized  search  algorithm  for  the  optional  grid  search  and  gradient  decent 
algorithm. 

Table  4.  Algorithm  complexity  for  Shack-Hartmann  and  SWAT  Sensors 


Sensor 

Shack-Hartmann 

SWAT 

Type  of 

Operation 

Complexity 
n  =  image 
length 

o„() 

Complexity 
n  =  image 
length 

OnO 

Additions 

2n2+n-2 

On(n2) 

2n-2 

On(n) 

Multiplications 

n 

On(n) 

n 

0„(n) 

Divisions 

i 

0„(1) 

i 

0„(1) 
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Table  5.  Algorithm  complexity  for  Maximum-Likelihood  Sensor  Phase  I 


Operation 

Phase  I 

Non-Optimized  Search 

Complexity 
n  =  image  length 

N  =  #  to  interpolate 

On()(*,+) 

OnO 

Complexity 
n  =  image  length 

N  =  #  to  interpolate 

o„()(*,+) 

OnO 

Logical  Shift 

(n-l)(N-l) 

On(n)* 

On(N) 

1 

On(l) 

ON(0) 

Additions 

(n-l)(N-l) 

On(n)* 

On(N) 

(n-l)((n/2)N+l)+(nN+2) 

On(n2)* 

On(N) 

Multiplications 

0 

on(0) 

ON(0) 

(n/2)((n/2)N+l) 

On(n2)* 

On(N) 

Natural  Logs 

(n-l)N+l 

On(n)* 

On(N) 

0 

o„(0) 

ON(0) 

Divisions 

0 

on(0) 

ON(0) 

0 

on(0) 

ON(0) 

Log-Likelihood 

Computations 

0 

on(0) 

ON(0) 

(n/2)N+l 

On(n)* 

On(N) 

Table  6.  Algorithm  complexity  for  Maximum-Likelihood  Sensor  Phases  II  and  III 


Where 


g’  =  1  if  including  the  optional  grid  search,  0  otherwise 
g  =  the  number  of  points  in  the  optional  search  grid,  or  one 


Operation 

Phase  II 

Phase  III 

Complexity 
n  =  image  length 
g  =  #  of  Grid  points 

N  =  #  to  interpolate 

o„0(*,+) 

OnO 

Complexity 
n  =  image  length 
g  =  #  of  Grid  points 

N  =  #  to  interpolate 

o„()(*,+) 

OnO 

Logical  Shift 

0 

o„(0) 

ON(0) 

1 

o„(i) 

ON(0) 

Additions 

3(n/2)g 

On(n) 

ON(0) 

Best:  (n+l)(2-2g’+log2((n/(2g))N)) 
Worst:  (n+l)(2-2g’+21og2((n/(2g))N)) 

On(nlog(n))+ 

0N(log(N)) 

Multiplications 

(n/2)g 

On(n) 

ON(0) 

Best:  (n/2)(2-2g’+log2(  (n/(2g))N)) 
Worst:  (n/2)(2-2g’+21og2((n/(2g))N)) 

On(nlog(n))+ 

0N(log(N)) 

Natural  Logs 

0 

On(0) 

ON(0) 

0 

O„(0) 

ON(0) 

Divisions 

0 

o„(0) 

oN(0) 

0 

On(0) 

ON(0) 

Log-Likelihood 

Computations 

G 

On(l) 

ON(0) 

Best:  2-2g’+log2((n/(2g))N) 

Worst:  2-2g’+21og2((n/(2g))N) 

On(log(n))+ 

0N(log(N)) 
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Image  Size  (one  Side)  {pixels} 

Figure  2 1 .  Illustration  of  Log-Linear  Nature  of  ML  Algorithm  for  LGS 
To  illustrate  these  results  for  varying  image  sizes,  Figure  21  provides  MATLAB 
simulation  results,  using  floating-point  computations  on  an  Athlon  2600+  processor,  of 
average  time  required  for  Shack-FIartmann  sensor,  cent  on  the  plot,  SWAT  wave-front 
sensor,  and  the  optimized  maximum-likelihood  algorithm  (mliw  with  linear  and 
MATLAB ’s  cubic-spline  interpolation)  for  different  image  sizes.  The  pre-compute  phase 
requires  the  same  computational  complexity  for  that  of  SWAT  estimation  for  additions 
but  subsequently  includes  the  additional  time  for  computing  logarithms  as  well.  The 
non-optimized  log-likelihood  search  algorithm  clearly  exceeds  the  computational 
complexity  of  both  centroiding  sensors;  however,  the  optimized  algorithm  is  only  slightly 
more  complex  compared  to  the  SWAT  wave-front  sensor.  Provided  the  pre-compute 
time  occurs  off-line  or  over  the  span  of  a  few  estimates,  the  optimized  algorithm 
performs  well,  with  0(n  log(n))  complexity,  which  is  far  superior  to  the  non-optimized 
search  algorithm. 
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4.2.5.  Possible  Improvements 

Through  investigation  of  the  requirements  and  assumptions  for  the  optimized 
search  algorithm  to  perform  properly,  a  few  improvements  are  possible  given  thorough 
understanding  of  the  algorithm  and  the  data  inputs  required.  This  optimized  algorithm 
assumes  that  the  interpolated  and  projected  true  image,  the  log  of  this  projection,  and  the 
observed  data  are  available  in  memory  for  multiple  accesses,  with  any  pre-computation 
completed  before  estimation  begins.  Overall  performance  improvement  may  be  possible 
by  computing  on-the-fly  interpolation  of  the  true  image,  which  also  requires  on-the-fly 
logarithms,  allowing  for  a  slightly  decreased  external  storage  space,  greatly  decreased 
data  access  times,  but  more  computations  overall.  If  implemented,  on-the-fly 
interpolation  could  lead  to  automatic  interpolation  level,  or  sub-pixel  step  size,  selection 
that  could  further  reduce  the  number  of  computations  required  for  the  search  algorithm. 

Although  this  optimized  search  algorithm  is  computationally  efficient  for  the 
number  of  log-likelihood  calculations  performed,  any  other  search  or  sort  algorithm 
capable  of  exploiting  new  hardware  technologies  and  possibly  breaking  apart  the  log- 
likelihood  calculation  itself  may  improve  upon  this  search  algorithm.  For  much  larger 
image  sizes,  and  possibly  joint  estimation  in  two  dimensions,  an  artificial  intelligence 
algorithm  could  provide  faster  results. 

4.2.6.  Limitations 

The  pre-computation  of  the  interpolation  and  logarithm  of  the  true  image  could 
significantly  hinder  the  effectiveness  of  this  wave-front  sensing  algorithm  depending  on 
the  nature  of  the  observed  image,  as  this  true  image  requires  updates  to  prevent  changes 
in  light-level,  orientation,  and  contrast  from  affecting  the  outcome  of  the  log-likelihood. 
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Additionally,  the  implemented  version  of  the  algorithm  requires  the  image  size  be  a 
multiple  of  four,  which  reduces  the  possible  array  sizes  for  a  given  AO  hardware  setup. 

4.3.  Implementation  Strategy 

4.3.1.  System  Layout 

The  benefit  of  a  maximum-likelihood  vector-correlating  wave-front  sensor  is  the 
ability  to  perform  tracking  and  wave-front  sensing  for  AO  using  existing  technologies 
with  hardware  requirements  that  typical  systems  already  meet.  As  long  as  the  imaging 
device  meets  or  exceeds  Nyquist  sampling  criteria,  and  the  image  provided  contains 
twice  or  greater  the  number  of  pixels  needed  for  the  search  window  with  appropriate 
light-level  and  contrast,  nearly  any  imaging  system  can  use  this  algorithm  as  a  tracking  or 
wave-front  sensor.  Since  an  appropriate  hardware  layout  exists  (lens,  array,  CCD,  etc) 
for  this  modern  wave-front  sensor,  this  section  focuses  on  the  efficiency  of  the  charge- 
coupled  device  (CCD)  array,  the  setup  of  the  processing  unit  hardware,  and  method  to 
implement  the  maximum-likelihood  estimation  algorithm  [5], 

The  size  of  the  image  follows  from  the  shift  detection  requirements;  however, 
exceeding  this  can  greatly  reduce  the  bias  and  overall  error  in  estimates.  Thus  an  image 
larger  than  thirty-two  by  thirty-two  pixels  is  the  recommended  size  for  tracking  and 
wave-front  sensing  of  plus  or  minus  four-waves  of  tilt,  corresponding  to  a  CCD  array  of 
160  or  greater  pixels  for  an  array  of  five-by-five  wave-front  sensors.  An  important 
possibility  to  consider  if  timing  constraints  pennit  is  multiplexing  the  x  and  y  dimension 
estimates  in  time,  allowing  for  twice  the  light  per  dimension  and  requiring  only  once 
CCD  array  instead  of  two.  If  the  CCD  array  has  a  fast  analog  to  digital  converter,  it 
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could  interpolate  the  data  directly  for  searching  of  shifts  in  the  observed  image  instead  of 
the  true  image,  which  offloads  this  requirement  from  the  pre-compute  phase. 

The  basic  setup  remains  the  same  except  for  the  above-mentioned  improvements 
and  the  specification  of  a  field  programmable  gate  array  (FPGA)  programmed  to  execute 
the  pre-compute  phase  and  search  algorithm.  The  targeted  architecture  for  this 
application  is  an  Altera  Cyclone  II  FPGA,  which  is  not  the  fastest  or  most  complex 
FPGA,  but  it  includes  some  basic  modules  to  aid  in  computation,  a  reasonable  package 
size,  and  a  greatly  reduced  cost  [7].  The  development  environment  for  this  device  is  free 
and  allows  the  researcher  to  validate  implemented  algorithms  in  the  targeted  environment 
with  a  simple  and  easy  to  use  interface. 

4.3.2.  VHDL  Implementation 

The  approach  to  algorithm  implementation  in  this  Very  High  Speed  Integrated 
Circuit  (VHSIC)  Hardware  Description  Language  (VHDL)  follows  a  basic  strategy  to 
break  the  problem  into  smaller  parts  and  implement  them  one  at  a  time,  which  requires  an 
extensible  state-machine  design  with  robust  transitions  and  outputs.  Abstraction  is 
critical  to  reduce  repeated  computations  to  manageable  modules,  prevent  excessively 
long  code,  and  allow  for  easier  debugging.  The  above  algorithm  reduces  easily  to  the 
Moore  state-machine  in  Figure  22,  allowing  for  robust  transitions,  controlled  outputs,  and 
extensions  to  include  further  options.  In  the  figure  the  symbol  “!”  indicates  negation,  the 
states  with  a  label  followed  by  a  colon  occur  in  the  indicated  arrangement  and  the  states 
outlined  in  dashes  are  extensions  to  this  diagram  currently  operating  as  direct 
connections  to  the  next  labeled  state  in  the  diagram. 
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Figure  22.  Moore  State-Machine  for  Implementation 
The  implementation  uses  32  bit  integer  rather  than  floating-point  computations  to 
increase  speed  and  reduce  complexity;  therefore,  to  preserve  accuracy  of  results  the  pixel 
values  of  the  true  image  and  the  log  of  such  values  simply  scale  up  by  the  required 
accuracy,  typically  four  decimal  places  for  the  simulated  light-levels,  allowing  correct 
computation  of  the  log-likelihood.  The  final  output  is  also  an  integer  value,  which  is  a 
linear  index  into  the  up-sampled  true  image  vector  and  the  estimated  shift  computes 
directly  from  Equation  6 1 . 
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Shift  = 


(61) 


Result 

_ jV _ 

s 

Shift,  Estimated  Shift  Value  (pixels) 

Result,  VHDL  Module  Result  (array  Position) 

IvuL,  Up-Sampled  True  Image  Length  (photons) 

S,  Amount  Multiplied  by  Nyquist  for  Over-Sampling  (unitless) 

[-1],  Optional  Subtraction  for  Indexing  Starting  at  One 
The  PreCompute  and  Compute  states  both  provide  loop  control  for  searching  over  the 
log-likelihood  values,  and  currently  perfonn  tasks  in  a  clocked-combinatorial  fashion; 
however,  these  nodes  could  control  additional  state  machines  for  further  robust  operation. 
Not  shown  in  the  diagram  is  a  separate  module  designed  specifically  to  compute  the  log- 
likelihood  values  providing  its  own  loop  control  to  step  through  the  windowed  images 
and  perform  computations.  This  does  not  implement  the  interpolation  or  logarithmic 
portion  of  the  maximum-likelihood  optimized  algorithm;  however,  preliminary  fitting 
and  timing  results  indicate  that  the  entire  algorithm  requires  less  than  four  percent  of  the 
logic  blocks  available  in  the  FPGA,  while  operating  at  a  conservative  clock  frequency  of 
33  MHz.  This  design  performs  a  single  log-likelihood  computation  every  twenty  clock 
cycles,  which  drives  a  total  computation  time  of  less  than  13.4  ps  for  a  single  estimate  on 
a  32  pixel  image  allowing  for  multiplexing  of  a  single  estimation  algorithm  for  seventy- 
five  sensors  at  an  update  rate  of  1,000  Hz.  This  speed  analysis  indicates  more  than 
adequate  temporal  efficiency  and  basic  simulation  results  appear  in  the  next  Chapter, 
while  the  VHDL  code  implementing  this  phase  of  the  algorithm  is  in  Section  4.2.3. 
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4.3.3.  Extensions  to  Implementation 

Including  the  interpolation  and  logarithm  in  the  FPGA  would  simplify  the  system 
greatly  by  allowing  a  single  infonnation  bus  from  the  CCD  to  the  FPGA,  reducing  the 
number  of  pins  required  and  lowering  overhead.  Additionally,  converting  the  loop 
control  statements  within  the  PreCompute  and  Compute  states  to  small  state  machines 
would  further  reduce  the  complexity  of  the  YHDL  code,  providing  extensibility  and 
versatility  to  the  system. 

4.4.  Summary 

This  Chapter  provides  the  analyses  for  three  major  investigations  in  tracking  and 
wave-front  sensing;  development  of  a  CRLB  for  any  tracking  or  wave-front  estimation 
techniques,  optimization  of  the  known  vector-projection  maximum-likelihood  algorithm, 
and  hardware  implementation.  To  validate  the  new  approaches  and  ensure  actual  results 
match  the  expectations  given  by  analysis  and  theory,  numerous  simulations  and 
parameterizations  in  the  next  chapter  cover  the  performance  of  operation  in  different 
environments. 
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V.  Results  and  Discussion 


This  Chapter  brings  together  the  previous  work  to  illustrate  in  an  intuitive  fashion 
the  perfonnance,  feasibility,  and  completeness  of  the  research  and  analysis  perfonned. 
First  presented  are  the  detailed  results  of  simulation  for  the  two-dimensional  Gaussian,  or 
LGS,  which  serve  to  narrow  the  range  of  interest  before  displaying  aggregate  results  for 
the  tracking  and  other  wave-front  sensing  applications.  Also  included  are  the  synthesis 
and  simulation  results  for  the  VHDL  implementation  of  the  maximum-likelihood 
optimized  algorithm.  As  indicated  in  the  Chapter  on  modeling,  the  simulation  step  size, 
0.1  pixels,  is  larger  than  the  search  step  size  as  set  by  the  CRLB  unless  indicated  different 
for  a  particular  plot.  The  total  number  of  noisy  trials  selected  for  each  realization  is  100 
trials-per-shift-per-dimension,  with  the  statistics  for  each  dimension  averaged  together, 
unless  indicated  different  for  a  unique  data  set. 

5.1.  LGS  /  2D  Gaussian  -  Modeling  and  Simulation  Results 

The  laser  guide  star  is  a  unique  case  in  that  the  analysis  providing  the  CRLB  uses 
this  type  of  image  for  the  bound  on  variance  for  any  estimator,  and  not  only  confirms  the 
design  choice  based  on  the  CRLB  but  also  provides  evidence  that  the  bound  is  indeed  a 
lower  bound  on  variance  for  wave-front  sensing.  The  parameters  used  for  modeling  the 
Gaussian  are  a  light-level,  C,  of  300  photons  and  a  standard  deviation,  o,,  of  two  pixels  as 
this  is  the  approximate  width  of  a  diffraction  limited  PSF. 
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5.1.1.  Effects  of  Image  Proj  ection 

To  provide  estimates  for  two-dimensions,  projection  of  the  Gaussian  image,  and 
any  image,  requires  either  an  equally  divided  light-level  or  an  equally  divided  time 
interval,  both  of  which  can  result  in  a  reduced  light-level  and  reduced  SNR.  Figure  23 
and  Figure  24  illustrate  the  effects  of  dividing  the  light  between  two  sensors  using  the 
whole  image  for  the  Shack-Hartmann  wave-front  sensor,  or  cent  corresponding  to  the 
first  CRLB,  and  projected  images  at  half  intensity  for  the  SWAT  wave-front  sensor 
corresponding  to  the  second  CRLB.  For  small  window  sizes,  the  edge  effects  not 
accounted  for  in  the  derivative  for  the  CRLB  appear  by  the  sensors  dipping  beneath  this 
bound;  however,  the  larger  16  pixel  image  illustrates  this  is  a  true  lower  bound  for  zero- 
shift  estimation  and  any  shifts  within  approximately  six  pixels  of  either  window  edge. 
This  data  represents  performance  results  without  a  background,  Bg;  and  this  parameter 
greatly  affects  the  centroid  algorithm  as  presented  later  in  this  chapter  in  Sub-Section 
5.3.3. 


Mean  Square  Error  for  Centroiding  Estimators 
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Shift  Value  {pixels} 


Figure  23.  Variance  of  Centroiding 
Algorithms  and  CRLB 


Figure  24.  Mean  Square  Error  of 
Centroiding  Algorithms 
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Since  the  MSE  is  the  bias  squared  plus  the  variance  as  indicated  by  Sub-Section  3.5.3, 
the  CRLB  is  still  a  lower  bound  on  this  statistic  as  well;  therefore,  MSE  plots  also 
include  the  CRLB,  with  both  plots  computed  over  10,000  trails-per-shift  in  each 
dimension  for  this  simulation.  The  MSE  indicates  that  although  the  variance  is  half  for 
the  Shack-Hartmann  wave-front  sensor,  the  bias  remains  constant  regardless  of  light- 
level,  and  this  bias  overwhelms  the  noise  error,  making  the  estimation  perfonnance  of  the 
Shack-Hartmann  and  SWAT  sensors  equally  poor  for  such  a  small  image  size. 

Given  that  the  Shack-Hartmann  wave-front  sensor  is  not  a  projection-based 
estimation  technique,  it  is  not  feasible  to  search  over  plus-or-minus  four  waves  of  tilt  in  a 
reasonable  amount  of  time.  The  SWAT  wave-front  sensor  is  eight-times  faster  in  reading 
data  for  an  eight  pixel  image,  which  only  covers  plus-or-minus  3.5  pixels  or  1.75  waves 
of  tilt,  eliminating  the  typical  Shack-Hartmann  centroiding  wave-front  sensor  from 
further  discussion. 

5.1,2.  Detailed  Bias 

As  mentioned  previously,  the  bias  indicates  the  best-case  operation  of  a  sensor 
and  is  difficult  to  remove;  therefore,  a  minimal  bias  is  ideal  for  any  type  of  sensor. 

Figure  25  displays  the  SWAT  and  maximum-likelihood  wave-front  sensors  for  the 
minimum  required  image  sizes  to  search  plus-and-minus  four  waves  of  tilt.  This 
simulation  uses  a  step  size  of  0.001  pixels  to  capture  the  quantization  error  caused  by  the 
CRLB  set  search  step  size  of  approximately  0.008  pixels. 
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Figure  25.  Absolute  Bias  v  Shift  for  LGS 

The  most  important  conclusions  to  draw  stem  from  the  clear  improvement  in  bias  of  the 
maximum-likelihood  sensor,  which  windows  the  data  and  the  true  image  before  shifting 
the  image,  over  the  SWAT  wave-front  sensor  and  the  nearly  order  of  magnitude 
difference  between  the  linear,  subscript  1,  and  cubic-spline,  subscript  c,  interpolated 
algorithms. 

The  expected  quantization  noise  of  approximately  0.004  pixels  due  to  the  search 
algorithm’s  interpolation  step  size,  which  the  CRLB  set,  is  nearly  perfect  with  the  cubic- 
spline  interpolated  algorithm.  This  graph  indicates  the  unbiased  nature  of  the  maximum- 
likelihood  algorithm  for  shifts  in  the  window,  provided  the  search  step  size  is  arbitrarily 
small  while  making  use  of  an  accurate  interpolator,  further  illustrating  that  the  CRLB  is  a 
correct  lower  bound  for  this  model  in  noise. 
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5.1.3.  Detailed  MSE  and  VAR 


Also  indicated  previously,  the  variance  illustrates  the  performance  of  the 
algorithm  in  noise,  and  tuning  of  the  sensor  without  changing  the  light  statistics  cannot 
reduce  or  eliminate  this  error.  The  mean  square  error  combines  the  variance  with  the  bias 
for  an  overall  performance  picture  as  shown  in  Figure  27,  and  the  variance  in  Figure  26, 
both  of  which  have  the  CRLB  with  half  the  light  intensity  for  the  projected  image  sensors 
and  10,000  trails-per-shift  in  each  dimension. 

The  key  points  to  derive  from  these  plots  include  the  clear  affect  of  bias  shown  in 
the  MSE,  making  some  estimated  shifts  useless  for  the  SWAT  sensor,  and  the  apparent 
noise  rejection  efficiency  of  all  sensors  for  relatively  small  shifts.  The  linear 
interpolator’s  error  appears  as  a  small  ripple  in  both  the  MSE  and  the  VAR,  but  nearly 
matches  the  cubic-spline  interpolator  with  respect  to  noise  rejection. 


Figure  26.  Variance  v  Shift  for  LGS  Figure  27.  Mean  Absolute  Error  v  Shift  for 

LGS 
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Although  it  does  not  affect  a  wave-front  correcting  AO  system,  it  is  important  to  note  that 
the  maximum-likelihood  sensor  does  not  always  return  the  edge  of  the  window  for  shifts 
significantly  outside  of  the  window,  which  is  a  direct  result  of  more  than  half  of  the 
object  moving  out  of  the  field  of  view.  Since  MSE  combines  the  variance  and  bias  of  a 
sensor  for  a  complete  picture,  further  results  only  display  the  mean  square  error. 

5.1.4.  Image  Size  Analysis 

Image  size  detennines  the  size  of  the  window  and  the  possibility  of  the  object 
moving  out  of  the  field  of  view,  and  although  the  minimum  image  sizes  are  perhaps  good 
enough  for  wave-front  sensing,  a  search  over  larger  images  sizes,  as  in  Figure  28  using 
10,000  trials-per-shift-per-dimension,  reveals  better  performance  for  all  sensors.  The 
following  results  include  averages  over  the  previously  displayed  regions,  with  a  solid  line 
on  the  graph  indicating  an  average  over  plus-or-minus  one-half  wave  of  tilt  while  a  dot- 
dashed  line  indicates  an  average  over  plus-or-minus  four- waves  of  tilt  or  the  entire  size  of 
the  window,  whichever  is  smaller.  Included  for  reference  for  this  plot  only  is  the  non¬ 
windowing  model  that  requires  the  functional  form  of  the  true  image  but  allows  better 
performance  with  a  smaller  image  size,  serving  to  illustrate  that  a  regular  windowing 
model  approaches  this  performance  with  an  image  16  pixels  larger.  The  overall  statistics 
indicate  an  image  size  of  32  pixels  is  optimal  for  the  non-windowing  maximum- 
likelihood  sensor,  mlnw  with  perfect  interpolation  and  extrapolation  using  the  Gaussian 
form  as  indicated  by  a  subscript  p.  An  image  size  of  48  pixels  allows  the  windowing 
models  that  do  not  require  explicit  knowledge  of  the  image  outside  of  the  window  to 
function  with  the  same  performance. 
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Figure  28.  Mean  Square  Error  v  Image  Size  for  LGS 
These  image  sizes  are  the  recommended  values  provided  the  hardware  and  subsequent 
readout  speeds  are  complementary;  however,  simulation  results  present  the  minimum 
image  sizes  to  illustrate  the  differences  and  trends  in  the  sensor  models. 

5.1.5,  Sampling  Analysis 

Over-sampling  of  an  image  aids  an  interpolator  by  providing  actual  data  in- 
between  pixels;  however,  it  increases  the  size  of  the  CCD  array  and  decreases  the  amount 
of  light  per-pixel  effectively  reducing  the  SNR;  therefore,  this  is  beneficial  only  if  the 
decrease  in  bias  is  greater  than  the  potential  increase  in  variance.  The  average  statistics 
in  Figure  29  indicate  that  over-sampling  aids  the  linear  interpolator  somewhat  and  harms 
the  cubic-spline  interpolator  to  a  lesser  extent,  with  the  exception  of  the  maximum- 
likelihood  data  shifting,  mldw,  and  swat  sensors,  which  perform  relatively  worse. 
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Multiple  of  Nyquist  Sampling  {unitless} 


Figure  29.  Mean  Square  Error  v  Sampling  Rate  for  LGS 
Although  the  corresponding  decrease  in  SNR  changes  the  variance  only  slightly,  the 
MSE  indicates  a  nearly  constant  trend  in  error  with  the  bias  driving  the  results.  These 
results  also  indicate  the  best  setup  for  optimal  performance  of  these  sensors  is  sampling 
at  or  slightly  beyond  Nyquist,  as  the  changes  in  overall  error  are  relatively  small. 

5.1.6,  Background  Intensity  Analysis 

Independent  of  image  size;  there  can  be  stray  light  in  the  optics  setup  as  well  as 
stray  electrons  in  the  CCD  array  itself,  which  the  sensor  sees  as  a  background  in  the 
image  lowering  the  overall  contrast  and  decreasing  the  SNR.  The  following  results  in 
Figure  30  show  the  importance  of  minimizing  the  background  for  the  centroid-based 
sensor  for  estimation  beyond  one-half  wave  of  tilt,  and  the  poor  performance  of  searching 
by  interpolating  and  shifting  the  observed  data  for  the  maximum-likelihood  algorithm. 
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Figure  30.  Mean  Square  Error  v  Background  Intensity  for  LGS 
The  windowing  maximum-likelihood  sensor  that  estimates  the  shift  by  shifting  the  image 
has  a  good  tolerance  to  background,  although  the  growth  rate  of  the  MSE  for  this  sensor 
still  appears  to  be  a  slow  exponential,  indicating  higher  background  would  disrupt  the 
performance  of  this  sensor  significantly  as  well.  As  the  SWAT  and  mldw  sensors  do  not 
perform  as  well  as  the  typical  maximum-likelihood  sensor,  the  discussion  will  no  longer 
include  the  centroiding  sensor  and  simulations  will  not  include  the  mldw  sensor. 

5.1.7.  Light-Level  Analysis 

The  CRLB  indicates  that  an  efficient  sensor  could  perform  better  given  higher 
light-levels,  and  Figure  3 1  shows  the  improved  performance  relationship  for  the 
maximum-likelihood  sensor  for  higher  light-levels.  As  noise  levels  decrease  below  the 
bias  of  the  sensor,  the  bias  drives  the  remaining  error  in  the  model 
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Mean  Square  Error  v  Light-Level 


Figure  3 1 .  Mean  Square  Error  v  Light-Level  for  LGS 
The  SWAT  wave-front  sensor  noticeably  demonstrates  this  trend  over  an  average  of  four 
waves  of  tilt.  As  light  levels  increase,  eventually  all  models  would  suffer  the  same  bias- 
limiting  effect;  however,  the  maximum-likelihood  sensors  have  the  advantage  that  a  finer 
search  step  size  driven  by  the  CRLB  allows  them  to  remain  efficient  in  noise  for  higher 
light-levels. 

5.2.  Tracking  Extended  Object  -  Modeling  and  Simulation  Results 

Possibly  the  most  interesting  simulation  of  this  chapter  is  tracking  an  extended 
object  and  these  results  exercise  the  extended  object  model  as  described  in  Sub-Section 

3.2.2.  There  are  two  difficulties  with  extended  objects;  first,  the  object  typically  fills  the 
field  of  view,  lowering  the  overall  contrast  of  the  image  and  second,  as  the  object  shifts 
new  information  enters  the  field  of  view,  which  can  cause  a  sensor’s  performance  to  drop 
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and  bias  to  be  asymmetric.  The  contrast  ratio  of  the  original  Hubble  image  before  proper 
band-limiting  is  approximately  four-to-one;  therefore,  these  results  should  suffice  for  any 
similar  object  with  a  similar  telescope  OTF  and  noise  statistics.  The  performance  graphs 
include  only  the  maximum-likelihood  sensor  as  a  centroid  algorithm  cannot  cope  with 
new  information  entering  the  scene  or  an  object  with  multiple  peaks  in  the  image. 
Additionally,  the  dimensions  are  separate  to  illustrate  the  performance  for  a  low-contrast 
orientation  and  an  orientation  containing  the  full  contrast  of  the  object,  which  indicates 
that  proper  tuning  of  the  sensor  or  optics  must  occur  before  tracking  is  possible. 

5.2.1.  Image  Size  Analysis 

As  opposed  to  Figure  30  for  the  LGS  case,  Figure  32  illustrates  an  increase  in 
required  image  size  for  acceptable  tracking  performance. 


Figure  32.  Mean  Square  Error  v  Image  Size  for  Hubble  Tracking 


94 


This  graph  also  indicates  that  the  CRLB  derived  for  the  LGS  using  the  same  parameters 
is  indeed  a  lower  bound  for  tracking  with  this  particular  image  of  the  Hubble  telescope  as 
2g,  /  C  =  1x10'  -pixels  .  The  subsequent  simulation  results  show  only  the  recommended 
image  size  of  48-by-48  pixels,  except  for  the  sampling  analysis  to  show  a  complete 
sampling  range. 

5.2.2.  Sampling  Analysis 

Figure  33  indicates  that  over-sampling  does  not  affect  the  performance  of  the  sensor. 
Again,  the  recommendation  is  to  sample  at  or  slightly  beyond  Nyquist  sampling  criteria 
to  avoid  aliasing  of  the  true  image. 


Multiple  of  Nyquist  Sampling  {unitless} 


Figure  33.  Mean  Square  Error  v  Sampling  Rate  for  Hubble  Tracking 
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5.2.3.  Background  Intensity  Analysis 

As  an  extended  object  already  has  a  background,  addition  of  further  background 
light  only  reduces  the  contrast  and  produces  small  changes  overall  in  the  performance  of 
the  maximum-likelihood  sensors  using  either  type  of  interpolation.  Figure  34  illustrates 
the  reasonable  tolerance  to  background  in  the  y-dimension  and  the  same  tolerance,  but 
consistently  poor  performance  for  the  low-contrast  x-dimension. 

5.2.4.  Light-Level  /  Total  Intensity  Analysis 

Much  like  the  laser  guide  star,  the  maximum-likelihood  sensor  performs  better 
with  higher  light-levels  as  shown  in  Figure  35;  and  the  roughness  of  the  curve  here  is  due 
to  the  limited  number  of  trails  for  such  a  large  bias  and  variance. 
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Figure  34.  Mean  Square  Error  v  Background  Intensity  for  Flubble  Tracking 
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Figure  35.  Mean  Square  Error  v  Light  Level  for  Flubble  Tracking 


5.3.  WFS  Extended  Object  -  Modeling  and  Simulation  Results 

Adaptive  optics  systems  attempt  to  correct  atmospherically  induced  wave-front 
distortions  for  any  object  the  researcher  wishes  to  view;  therefore,  a  simulation  of  a 
wave-front  sensing  application  on  the  Hubble  provides  further  insight  into  the  utility  of 
this  sensor,  and  the  possibility  of  tracking  and  imaging  with  the  same  telescope  and 
sensor  setup. 

5.3.1.  Image  Size  Analysis 

Figure  36  again  confirms  that  wave-front  sensing  over  plus-or-minus  four  waves 
of  tilt  is  possible  at  the  minimum  image  size,  but  the  recommended  image  size  for  the 
maximum-likelihood  wave-front  sensor  remains  a  reasonable  48  pixels  in  length. 
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Figure  36.  Mean  Square  Error  v  Image  Size  for  Flubble  WFS 


The  maximum-likelihood  wave-front  sensors  maintain  good  performance  for  larger 
image  sizes,  while  the  centroiding  algorithm’s  performance  decreases  as  the  image  size 
increases  due  to  the  greater  incorporation  of  the  image’s  natural  background.  This  image 
also  indicates  that  the  CRLB  is  a  good  approximation,  given  the  parameter  a;  is  an 
approximation  for  the  average  image  shape  in  both  dimensions  of  the  image. 

5.3.2.  Sampling  Analysis 

This  simulation  results  in  Figure  37  indicate  the  same  increase  in  performance 
with  respect  to  over-sampling  for  the  linear  interpolated  sensor,  and  the  same  slight 
decrease  in  performance  for  the  centroiding  and  cubic-spline  interpolated  sensor. 
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Figure  37.  Mean  Square  Error  v  Sampling  Rate  for  Flubble  WFS 


5.3.3.  Background  Intensity  Analysis 

Despite  the  inherent  background  of  the  Hubble  image,  the  maximum-likelihood 
sensor  continues  to  demonstrate  strong  tolerance  to  further  background  light  in  Figure  38, 
with  very  reasonable  mean  squared  error  performance.  The  swat  wave-front  sensor 
degrades  at  a  slower  rate  than  seen  with  the  laser  guide  star  due  to  the  higher  light  level 
of  this  object. 

5.3.4.  Light-Level  /  Total  Intensity  Analysis 

As  the  CRLB  predicts,  Figure  39  indicates  all  sensors  benefit  from  the  increased 
light-levels;  however,  the  increased  light-level  appears  to  aid  the  cubic-spline  interpolator 
slightly  more  than  the  linear  interpolator  due  to  the  linear  interpolator’s  added  error. 
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Figure  38.  Mean  Square  Error  v  Background  Intensity  for  Flubble  WFS 


Figure  39.  Mean  Square  Error  v  Fight-Fevel  for  Hubble  WFS 
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These  simulations  illustrate  excellent  performance  for  the  vector-projection  maximum- 
likelihood  wave-front  sensor  as  well  as  validate  the  sensor  is  efficient  with  respect  to  the 
CRLB  for  larger  image  sizes.  The  CRLB  also  holds  for  other  image  types,  given  an 
appropriate  estimate  of  a,. 

5.4.  Implementation  /  VHDL  Simulation  Results 

The  simulation  results  for  YHDL  represent  the  estimated  speed  of  operation  for 
the  hardware  implementation  of  the  sensor;  however,  these  results  do  not  fully  test  the 
implementation  and  further  validation  should  proceed  before  reliance  on  this  realization 
of  the  sensor. 

5.4.1.  Targeted  Device  Resource  Summary 

The  following  results  in  Table  7  and  Table  8  indicate  the  usage  of  the  logic 
elements  in  the  FPGA  and  the  maximum  clock  speed  that  the  device  could  operate  using 
the  Cyclone  II  EP2C70F89618  processor. 


Table  7.  Resources  Required  for  Synthesis 


Resource 

Number  Used 

Number  Available 

Percent  Used 

Logic  Elements 

2,447 

68,416 

<4% 

Registers 

491 

I/O  Pins 

228 

622 

37% 

Memory  Bits 

0 

1,152,000 

0% 

Embedded  9  bit 

Multipliers 

6 

300 

2% 
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Table  8.  Preliminary  Timing  Analysis 


Timing  Type 

Value 

Clock 

53.79  MHz  (18.59  ns) 

Worst-Case  Setup  Time  (tsu) 

23.381  ns 

Worst-Case  Hold  Time  (th) 

-1.225  ns 

Worst-Case  Clock  to  Output  (tco) 

14.21  ns 

Additionally,  the  preliminary  estimated  power  consumption  is  a  low  260  mW  allowing 
for  reduced  cooling  requirements  in  an  embedded  environment.  These  results  indicate 
that  a  clock  period  of  30  ns  is  a  safe  value  for  setup  and  hold  times  while  maintaining  the 
performance  of  the  estimation  algorithm,  which  results  in  a  total  compute  time  of  less 
than  13.4  ps,  allowing  for  time-division  multiplexing  of  one  module  for  up  to  75  different 
wave-front  sensors  at  an  update  rate  of  1,000  Hz. 

5.4.2,  Test  Bench  Simulation  Results 

This  simulation  serves  to  indicate  the  total  time  required  to  complete  one  log- 
likelihood  search  for  a  32-by-32  pixel  image,  and  does  not  indicate  complete  accuracy  or 
validate  the  implementation  outside  of  the  simple,  non-realistic  test  case  presented,  as  it 
has  only  44.8  %  test  coverage.  Figure  40  indicates  the  form  of  an  individual  log- 
likelihood  computation  and  is  a  zoomed  view  of  the  beginning  of  Figure  41,  which  shows 
the  overall  computation  for  a  nearly  worst  case  of  23  log-likelihood  computations. 
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Figure  40.  Zoomed  in  View  of  First  Portion  of  VFIDL  Simulation 
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Figure  41.  Full  View  of  VFIDL  Simulation 


These  results  indicate  the  total  required  operating  time  as  well  as  proper  operation 
for  the  simplistic  test  case  presented.  Complete  algorithmic  flow  testing  and  addressing 
verification  are  follow-on  research  areas  in  a  complete  implementation. 
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5.5.  Summary 

The  three  main  areas  of  research  results  presented  in  this  chapter  conclude  nearly 
half  of  the  developmental  life  cycle  of  a  project  including  theoretical  research, 
algorithmic  development,  and  hardware  implementation,  with  the  remaining  portion  of 
the  life  cycle  including  at  a  minimum  hardware  realization,  testing,  and  maintenance.  As 
an  additional  benefit,  further  simulation  results  parameterizing  the  sensors  over  a  wider 
range  appear  in  Appendix  D.  The  analysis  for  the  CRLB  provides  implementation  and 
validation  benefits,  while  the  developed  log-likelihood  search  algorithm  and 
implementation  provide  a  solid  foundation  and  proof  of  concept  for  implementation. 
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VI.  Conclusions 


From  the  results  and  discussion,  it  is  clear  that  the  vector-projection  maximum- 
likelihood  wave-front  and  tracking  sensor  is  quite  versatile  and  feasible  to  implement. 
This  research  intended  to  and  succeeded  in  providing  definitive  results  for 
characterization  and  implementation  of  the  maximum-likelihood  sensor  through  the  use 
of  applied  theory,  robust  modeling,  and  sound  implementation  techniques.  Combined 
with  the  power  of  a  fast  search  algorithm,  maximum-likelihood  could  become  the  new 
standard  in  embedded  estimation  techniques. 

6.1.  Key  Contributions 

Not  only  did  this  research  develop  a  Cramer-Rao  lower  bound  for  any  wave-front 
sensing  or  tracking  application  but  it  further  reduced  this  complex  theoretical  model  to  a 
simple  point  solution  providing  an  easy  method  to  predict  and  validate  real-life  results. 
The  results  of  the  bias  analysis  provides  targeted  implementation  information  as  well  as 
the  noise  statistics  which  further  solidified  the  utility  of  this  sensor  by  indicating  the 
greater  efficiency  compared  to  the  centroiding  algorithms  currently  in  use.  The 
optimized  search  algorithm  with  noise  rejection  capability  designed  for  a  concave-down 
log-likelihood  is  useful  for  any  maximum-likelihood  application  with  a  similar 
likelihood,  log-likelihood,  or  other  type  of  curve  to  seek  a  maximum  or  minimum  over. 
Finally,  the  development  method  and  implementation  of  this  optimized  algorithm  in 
hardware  reduces  the  life-cycle  time  greatly  for  current  use,  while  providing  a  road  map 
for  future  implementations  of  other  complex  search  algorithms  or  embedded  software. 
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6.2.  Lessons  Learned 


Several  difficult  areas  of  this  research  could  have  been  simpler  if  more  time  and 
care  had  focused  on  proper  modeling  of  the  environment  for  statistical  analysis  and 
simulation.  It  is  impossible  to  determine  if  an  algorithm  or  implementation  is  correct  if 
the  information  sent  to  the  element  of  research  is  not  correct,  as  only  unexpected  results 
occur.  In  addition,  reliance  on  modern  simulation  tools  often  led  to  further  problems  as 
the  implementation  of  these  tools  is  not  always  clear  and  may  provide  inaccurate  results 
in  limited  circumstances.  Outside  of  these  areas  of  difficulty,  the  research  went  quite 
smoothly  by  capitalizing  on  the  background  and  understanding  provided  by  a  thorough 
education. 

6.3.  Further  Research 

As  with  any  research  there  are  many  areas  in  which  improvements,  extensions, 
and  developments  are  available;  therefore,  this  section  concludes  with  the  possible 
follow-on  work  associated  with  the  areas  of  research  investigated. 

The  Cramer-Rao  lower  bound  only  includes  the  intensity  of  light  as  a  parameter, 
ignoring  contrast,  which  is  the  single  largest  contributor  to  reduced  perfonnance  for  an 
extended  object.  The  realization  of  an  inclusion  or  relation  to  contrast  would  provide 
even  greater  utility  to  the  simplistic  yet  effective  zero-shift  minimum  variance 
calculation.  Additionally,  the  CCD  estimates  each  pixel  in  an  image  during  the  capture 
process;  therefore,  a  more  accurate  bound  would  jointly  estimate  the  current  parameters 
with  each  pixel  intensity. 

The  vector-projection  maximum-likelihood  tracking  and  wave-front  sensor  has 
numerous  areas  of  investigation,  only  some  of  which  stem  directly  from  this  research.  It 
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may  be  possible  to  compute  a  closed-form  derivative  for  any  image,  allowing  for  single- 
step  computation  for  ML  estimation.  With  the  ability  to  dynamically  search  the  log- 
likelihood,  it  would  be  very  possible  to  perform  on-the-fly  interpolation  and  possibly 
even  automatic  interpolation  level  adjustment  based  on  the  results  of  the  optimized 
search  algorithm.  Another  important  area  of  investigation  is  the  use  of  this  search 
algorithm  for  joint  estimation  in  two-dimensions,  either  using  vector-projection  or  the 
entire  image,  or  other  techniques  such  as  phase  diversity  estimation.  An  enhancement  to 
the  tracking  application  would  be  the  inclusion  of  automatic  light-level  normalization 
between  estimations,  which  would  provide  greater  accuracy  and  better  performance. 

The  embedded  implementation  should  match  the  performance  of  a  simulation 
algorithm;  however,  further  state-machine  analysis  would  simplify  the  hardware  layout 
and  provide  a  more  robust  solution.  Further  testing  and  development  of  this  algorithm  to 
validate  the  design  could  allow  for  immediate  implementation.  Finally,  implementation 
of  the  required  interpolation  and  logarithmic  functions  would  increase  the  overall 
productivity  of  an  embedded  system  and  simplify  the  data  transfer  requirements  between 
different  systems  greatly. 

Even  without  the  above  extensions,  maximum-likelihood  has  a  promising  future 
in  the  arenas  of  tracking  and  wave-front  sensing  applications. 
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Appendix  A:  Mathematica  Verification  of  CRLB 


A.f.  Setup  of  Formulas 


ln[l ]:=  Off  [General: :  spell] 

c 

In[2]:=  i  [p_,  o_,  C_,  r_]  :=  —==  e  2<l2 

V  2  7C  v2 
«< 

lnp]:=  JLL[£_,  o_,  C_,  d_]  :  =  ^  (dLo»[i[/S,  o,  C,  r]]  -Log[d!]  -i[/S,  a,  C,  r]) 


In [4] :=  fpJl[p_,  a_,  C_,  d_]  :=  0£  JLL[fi,  a,  C,  d] 

fpJ2[P_,  a_,  C_,  d_]  :=  0CT  JLL[/5,  a,  C,  d] 

fpJ3[p_,  o_,  C_,  d_]  :=  St  JLL[p,  a,  C,  d] 


ln[7]:=  spJll[0_,  or_,  C_,  d_]  :=  dp  fpJl[0,  a,  C,  d] 

spJ12[p_,  a_,  C_,  d_]  :=  0c  fpJl[p,  a,  C,  d] 

spJ13[p_,  <r_,  C_,  d_]  :=  0t  fpJl[p,  a,  C,  d] 

spJ21[p_,  a_,  C_,  d_]  :=  dp  fpJ2[/3,  a,  C,  d] 

spJ22[p_,  <r_,  C_,  d_]  :=  0,  fpJ2[p,  a,  C,  d] 

spJ23[p_,  a_,  C_,  d_]  :=  ^  fpJ2[/i,  a,  C,  d] 

spJ31[/5_,  a_,  C_,  d_]  :=  0jq  fj>J3[p,  o,  C,  d] 

spJ32[p_,  a_,  C_,  d_]  :=  0ff  fpJ3[/J,  a,  C,  d] 

spJ33[p_,  a_,  C_,  d_]  :=  0fc  fi>J3[/5,  a,  C,  d] 


In[l6]:=  Jll[^_,  a_,  C_]  :=  -spJll[p,  a,  C,  d[r]]  /.  d[r]  -»i[p,  a,  C,  r] 

J12[0_,  a_,  C_]  :=  -spJ12  [f>,  a,  C,  d[r]]  /.  d[r]  ->i[P,  a,  C,  r] 

J13[/5_,  a_,  C _]  :=  -spJ13[p,  a,  C,  d[r]]  /.  d[r]  -»i[p,  ff,  C,  r] 

J21[p_,  a_,  C _]  :=  -spJ21[p,  a,  C,  d[r]]  /.  d[r]  ->![£,  a,  C,  r] 

J22[p_,  o_,  C  ]  :=  -spJ22[p,  a,  C,  d[r]]  /.  d[r]  -»i[p,  o,  C,  r] 

J23[0_,  a_,  C _]  :=  -spJ23[P,  a,  C,  d[r]]  /.  d[r]  ->i[p,  a,  C,  r] 

J31[0_,  o_,  C_]  :=  -spJ31[p,  a,  C,  d[r]]  /.  d[r]  ->i[/3,  a,  C,  r] 

J32[p_,  a_,  C_]  :=  -spJ32[p,  a,  C,  d[r]]  /.  d[r]  ->i[p,  a,  C,  r] 

J33[0_,  a_,  C _]  :=  -spJ33[0,  a,  C,  d[r]]  /.  d[r]  -»i[P,  a,  C,  r] 


In [25] :=  FP[P_,  a_,  C_,  d_] 


fpJl[P,  o,  C,  d] 
fpJ2[£,  a,  C,  d] 
fpJ3[P,  a,  C,  d] 


In [26]:=  SP[P_,  a_,  C_,  d_] 


spJll[P,  a,  C,  d]  s|)J12[p, 
spJ21[P,  a,  C,  d]  spJ22[]S, 
spJ31[p,  a,  C,  d]  spJ32[jG, 


,  C,  d]  spJ13[£,  a,  C,  d] 
,  C,  d]  spJ23[]8,  a,  C,  d] 
,  C,  d]  spJ33[£,  a,  C,  d] 


In [27]:=  FI[/5_,  a  ,  C_] 


Jll[p,  a,  C]  J12 [p,  a,  C]  J13[P,  a,  C] 

J21[0,  a,  C]  J22 Ip,  a,  C]  J23[0,  a,  C] 

J31[0,  a,  C]  J32[p,  <r,  C]  J33[p,  o,  C] 


In [28]:=  CRLB[p_,  o_,  C_]  :=  Inver se  [FI  [p ,  a,  C]] 
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A.2.  Computation  of  CRLB 


in [29]:=  Matri^orm[FiillSini>lify  [FP  [£,  a,  C,  d[r]]]] 

Out[29]//totatrixForm= 


(r^) 


(T-fi)* 

t  a1 


JT 


+t  d[r] 


£<j2 


(r-B)2 

<r-B)2 

«  2,2  (r2-2rfi*fi2-»2) 

a/T  C-£  e  t  at  ^  x  a/  (ft  d[r] 

<t5  +[** 


(* Verified*) 


In [30]:=  MatrixForjn[FullSiiTOlify[SP[/i,  a,  C,  d[r]]]] 


(* Verified*) 


Out  [30]//tofatrix  F  orm= 

Z_ 

Z„ 

2L 


(r-jB)2 

(r^S)2 

e  2  <'2  <r^) 

a/  t  C  r2-£->/  : 

Ce  2»!  (r^B)2  |  te  tv*  d[r] 

V’  . 

a/7*  <r*+j-<rt  a/  tn  a^A^at  ** 

Z - Jr=-*> 

Cr^B)2 

e  2  "2  (nB) 

(r^B)2  I 

Z2^  Ct2-t^/7  CrAv</7  CB2-3a/7  C  l2+4e  2"2  4/7  (ir2)3/2a[r] 

Cr^B)2 
e  £  <r2 

v 

-if*  t  r4VZ+4^/7 

tA^H.  <t5a/"^ 

Zl_Jr=-c* 

(isP)* 

e  tv*  ( -r+fi) 

aJT^  {'tyn 


_(£^) 1 

e  * (-r*+£  rjQ-jfi*+<T2) 


In [31]:=  MatrixF orm[FullSiui>lify[FI  a,  C]]] 


(*Verified*) 


Out  [3 1  ]//tofatrix  F  orm= 

-Z„ 

ZL 

ZL 


-ll=22! 

C  e  » »*  (.i-sf 

«]  Zn  vi'fvZ 

te  2  1,2  (  r^B)  (  -r2+*  r  jB^B2  V  ) 

,\J  <rS  a/72^ 

0-8)  2 

e  2  »2  (-r-t8) 

,/77  (»2)3'2 


(r^B)2 

I  e  Z  *2  Cr^Q)  (-r2+g  rj6^62+.;2) 

-v/77  T5'\/Z 

Ce  21,8  (»i-irM1-.i)i 

^  £  *  a6  a/"^ 

e  *  **  (-r*+£  rjQ-02+<X2) 

a>I  in  a^^fat 


_<r-^ 
e  *  v*  (-r-tJfi) 

a/"**  (**)*/* 

(r-g)2 

e  *  g*  (-r*+£  rj6-jfi2+«r2) 

a/  £  n  a 

(r-£)2 
e  £  IT2 

C  a/T*  a/"*2^ 


!tatrijfForm[FullSinf>lify[CRLB[£,  a,  C]]] 


(  *Incomi)  rehensUtle  *  ) 
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A.3.  Simplification  Setup  and  Solution 


In [32]:=  sJll[£_,  <7_,  C_]  :  =  assuming[o  >  0  44  r  e  Integers, 


sJ12[/S_,  «r_,  C_]  :  =  Sssuming[a  >  0  44  t  e  Integers, 


sJ13[/i_,  o_,  C_]  :  =  assuming [a  >  0  44  r  e  Integers, 


sJ21[/l_,  o_,  C_]  :  =  assuming [a  >  0  44  r  e  Integers, 


sJ22[p_,  o_,  C_]  :  =  assuming [o  >  0  44  r  e  Integers, 


sJ23[p_,  o_,  C_]  :  =  assumng[a  >  0  44  r  e  Integers, 


sJ31[£_,  o_,  C_]  :  =  assuming [o  >  0  44  r  e  Integers, 


sJ32[/i_,  o_,  C_]  :=  assuming [o  ■»  0  44  r  c  Integers, 


sJ33[p_,  a_,  C_]  :  =  assuming[o  >  0 44 r  e  Integers, 


I 

I 

I 

I 

I 

I 

I 

I 

I 


(.Z-P)2 

Ce  2^  <r-^fa] 

V  2  71  O4  V®2 

<z-0)2 

Cl  2V2  (r -/S)  (r2  -  2r/s+p2  -  o2) 
'/  2  71  O5  V®2 


(r-£)2 

e  2  <r-0) 

— ; - - dir] 

■'./  2  7i  (a2)3/2 
(t-P)2 

Ce  Ji2  (r  -/S)  (r2  -  2r  p+f?  -  o2) 
V  2  7T  O5  V®2 


(c-P)2 

Ce  202  (r2  -  2rp  +  p2-o2)2  ^ 

- - - - dr 

V  2  71  o6  V"®2 


(r-P)2 

e  si2  (r2  -  2r  p  +  p2  -  o2) 
V  2  7r  a3  Vo2 


dir] 


(C-P)2 

*  ^  <r-»«r] 

V  2  7T  (02)3/2 


(=A8)2 

e  2  d2  (r2  -  2  r  p  +  p2  -  a2) 
V  2  71  a3  Va2 


dir] 


O^)2 

e  2d2 

— - - dir] 

C  V  2  71  V"®2 


dir] 


dir] 


(sJll[0,  a,  C]  sJ12[p,  a,  C]  sJ13[0,  a,  C]  i 

sJ21[p,  o,  C]  sJ22[p,  o,  C]  sJ23[p,  o,  C] 

sJ31[/S,  a,  C]  sJ32[0,  a,  C]  sJ33[0,  a,  C]  J 

In [42]:=  sCRLB[/i_,  a,  C_]  :=  Inverse  [sFI  [p,  a,  C]] 

In [43]:=  MatrixForm[FullSini>lify[sFI  [£,  o,  C]]] 

Out  [43]PTiyfatrix  Form= 

f4  0  o 

<T* 

0  ^  0 

0  0  i 

C 

In [44]:=  MatrixForiTi[FullSiiTiJlifY[sCl*LB [0,  a,  C]]] 

Out  [44]//Tubtrix  F  orm= 

f-  0  0] 

C 

0  0 
t  c 

i0  0  Cl 
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Appendix  B:  MATLAB  Version  of  ML  Algorithm 


%  Maximum  Likelihood  Estimate  Using  Optional  Grid  &  Gradient  Decent 
Search 

%  function  a  =  mliwb  (Dv, Ivi, llvi, N, aMax, Grd) 

%  Dv  -  Data  Vector 

%  Ivi  -  Image  Vector  (Interpolated) 

%  llvi  -  Log  of  Interpolated  Image  Vector 

%  N  -  Interpolation  Points  Used 

%  aMax  -  Maximum  Search  Area  =  +/-  aMax 

%  Grd  -  Step  Size  for  Grid  Search  0  <  Grd  <  length (Dv)  to  Perform 

function  a  =  mliwb  (Dv, Ivi, llvi, N, aMax, Grd) 


Grd  =  f loor (Grd*N) /N; 
step  =  1/N; 


%  Ensure  Grid  Conforms  to  vecs 
%  Interpolations  Step  Size 


DvL  =  length (Dv); 

DvuL  =  (DvL-1 ) *N+1 ; 

IvuL  =  length (Ivi); 

IvL  =  (IvuL-1) /N+l; 
if  IvL  <  DvL 
a  =  NaN; 
return; 

elseif  IvL  <=  2*DvL 

sLim  =  (IvL/2-l)/2; 

else 


%  Data  Vector  Length 
%  Upsampled  Data  Vector  Length 
%  Interpolated  Image  Vector  Len 
%  Image  Vector  Length 
%  Check  for  incorrect  input  vectors 


%  Get  Maximum  Window  Size 


sLim  =  (DvL-1) /2; 

end; 

wLim  =  sLim  -  (DvL  -  2*sLim+l  -  1 ) / 2 ;  %  Get  Minimum  Window  Size 

if  aMax  >  sLim  %  check  aMax  is  Less  than  Maximum  Window 

aMax  =  sLim; 

SI  =  2*sLim  -  aMax; 

elseif  aMax  <  wLim  %  Check  aMax  is  Greater  than  Minimum  Window 

SI  =  sLim; 


else 


SI  =  2*sLim  -  aMax; 

end; 

iSg=  ( IvL-1 ) / 2 ; 

iSs=  (iSg-Sl+1-1) *N+1; 

iSe=  (iSg+Sl+1-1) *N+1; 

dSg=  (DvL-1 ) / 2; 

dSs=  dSg-Sl+1; 

dSe=  dSg+Sl+1; 

sSs=  (iSg-aMax+1-1) *N+1; 

sSe=  (iSg+aMax+1-1) *N+1; 


%  start  index  for  smaller  image 
%  end  index  for  smaller  image 

%  start  index  for  smaller  data 
%  end  index  for  smaller  data 

%  start  index  for  minimum  Shift  (in  window) 
%  end  index  for  maximum  Shift  (in  window) 


if  Grd  ==  0  | |  Grd*N  >  sSe-sSs  %  Determine  if  Grid  Search  is  Required 
noGrid  =  true; 

else 

noGrid  =  false; 

end; 


taL  =  zeros (1, IvuL) ; 


%  Allocate  Memory  (overkill) 
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if  noGrid  %  Just  Compute  Endpoints  (no  Grid) 

for  idx  =  [sSs  sSe]  %  Get  LL  Values  for  Endpoints 

taL(idx)  =  sum (Dv (dSs : dSe) . * . . . 

llvi ( (iSs:N:iSe)+ (IvuL-1) /2- (idx-1) ) . . . 

-Ivi ( (iSs :N: iSe)  +  (IvuL-1) /2- (idx-1) )  ,  2) ;  %  Shifts  Image 

end; 

if  taL(sSs)  <  taL(sSe)  %  Right-Side  is  Bigger 

max  =  sSe; 
sdx  =  2; 

else  %  Left-Side  is  Bigger 

max  =  s  S  s ; 
sdx  =  1; 

end; 

idx  =  sSs+f loor ( ( sSe-sSs ) /2 ) ;  %  Compute  Next  Search  Point 

else  %  Perform  Grid  Search 

init  =  sSs : Grd*N : sSe;  max  =  sSs; 

if  mod ( sSe-sSs , Grd*N)  ~=  0  %  Add  Last  Point  if  Necessary 

init  =  [init  sSe] ; 

end; 

for  idx  =  init  %  Quick  search  for  points  on  peak 

taL(idx)  =  sum (Dv (dSs : dSe)  . *  .  .  . 

llvi( (iSs:N:iSe)+ (IvuL-1) /2- (idx-1) ) . . . 

-Ivi ( (iSs:N:iSe)  +  (IvuL-1) /2- (idx-1) )  ,2)  ; 
if  taL(idx)  >  taL (max)  %  Replace  Maximum  if  Necessary 

max  =  idx; 

end; 

end; 

if  max  ==  sSs  %  If  Max  is  Left-Side 

taL(max+l)  =  sum (Dv (dSs : dSe) . * . . .  %  Get  Next  Value 

llvi ( (iSs:N:iSe) + ( IvuL-1 )  / 2- (max+1-1 ) )  . . . 

-Ivi ( (iSs :N: iSe) + (IvuL-1) /2- (max+1-1) ) , 2) ;  %  Shifts  Image 

if  taL (max)  <  taL(max+l)  %  Right-Side  Bigger 

sSs  =  max+1; 

if  max+Grd*N  <=  sSe  %  Decide  if  EndPoint 

sSe  =  max+Grd*N; 

else 

sSe  =  sSe; 

end; 

sdx  =  1; 

idx  =  sSs+f loor (( sSe-sSs ) /2 ) ;  %  Compute  Next  Search  Point 

else  %  Left-Side  is  Peak 

idx  =  -1;  %  Finish 

end; 

elseif  max  ==  sSe  %  If  Max  is  Right-Side 

taL(max-l)  =  sum (Dv (dSs : dSe) . * . . .  %  Get  Previous  Value 

llvi ( (iSs:N:iSe) + ( IvuL-1 ) / 2- (max- 1-1 ) ) . . . 

-Ivi (( iSs : N : iSe) +( IvuL-1 ) /2- (max-1-1 )), 2 ) ;  %  Shifts  Image 

if  taL (max)  <  taL (max-1)  %  Left-Side  Bigger 

if  mod ( sSe-sSs , Grd*N)  ==  0  %  Determine  Previous  Point 

sSs  =  max-Grd*N; 

else 

sSs  =  max-mod ( sSe-sSs , Grd*N)  ; 

end; 

sSe  =  max-1; 
sdx  =  2; 

idx  =  sSs  +  f loor (( sSe-sSs ) /2 )  ;  %  Compute  Next  Search  Point 

else  %  Right-Side  is  Peak 
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idx 


-1; 


Finish 


end; 

else  %  Max  is  in  Middle 

taL(max+l)  =  sum (Dv (dSs : dSe) . * . . .  %  Get  Next  Value 

llvi  (  (iSs:N:iSe)  +  (IvuL-1) / 2- (max+1-1) )  .  .  . 

-Ivi ( (iSs :N: iSe) + (IvuL-1) /2- (max+1-1) ) , 2) ;  %  Shifts  Image 

if  taL (max)  <  taL(max+l)  %  Right-Side  is  Bigger 

sSs  =  max+1; 

if  max+Grd*N  <=  sSe  %  Decide  if  EndPoint 

sSe  =  max+Grd*N; 

else 


sSe  =  sSe; 

end; 

sdx  =  1; 

idx  =  sSs  +  f loor (( sSe-sSs ) /2 )  ; 
else  % 

sSs  =  max-Grd*N; 
sSe  =  max; 
sdx  =  2; 

idx  =  sSs+f loor (( sSe-sSs ) /2 ) ; 

end; 

end; 

end; 


%  Compute  Next  Search  Point 
Left  Side  is  Bigger 


%  Compute  Next  Search  Point 


while  idx  >  0 

if  idx  <=  sSs 
idx  =  -1; 

else 


%  Check  for  Complete  Condition 
|  idx  >=  sSe  %  Check  to  See  if  we're  Done 

%  Finish 


taL ( idx) =sum (Dv (dSs : dSe) . * . . .  %  Compute  Next  LL  Value 

llvi  (  (iSs:N:iSe)  +  ( IvuL-1 ) /2- ( idx-1 ) )  .  .  . 

-Ivi ( ( iSs : N : iSe) + (IvuL-1) /2- ( idx-1) ) , 2 ) ;  %  Shift  Image 


if  taL ( idx)  > 

taL (max) 

%  New  Value  Bigger  than  Old  Max 

if  idx+1 

<  sSe 

%  Get  Next  Value 

taL ( idx+1 ) =sum (Dv (dSs : dSe) .  *  .  .  . 

llvi ( (iSs :N : iSe) 

+  (IvuL-1) /2-((idx-l)+l))  ... 

- 

Ivi ( (iSs :N : iSe) 

+ (IvuL-1) /2-( (idx-1) +1) ) ,2) ; 

end; 

dir  =  taL(idx)  <  taL(idx+l);  %  Store  Sign  of  Slope 

if  dir 

%  Store  Correct  Index 

idx  = 

idx+1 ; 

end; 

max  =  idx 

r 

%  Store  New  Max 

if  ~dir 

%  Shutter  Right  (Peak  Left) 

sSs 

=  s  S  s  ; 

sSe  = 

idx; 

sdx  = 

2; 

else 

%  Shutter  Left  (Peak  Right) 

sSs  = 

idx; 

sSe 

=  sSe; 

sdx  = 

1; 

end; 

else 

%  Slope  is  Toward  Current  Max 

max  =  max; 

if  sdx  == 

1 

%  Shutter  Right  (Peak  Left) 

sSs 

=  s  S  s  ; 

sSe  = 

idx; 

sdx 

=  1; 
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o\°  o\° 


else 


Shutter  Left  (Peak  Right) 


sSs  =  idx; 
sSe  =  sSe; 
sdx  =  2; 

end; 

end; 

idx  =  sSs+f loor ( ( sSe-sSs ) /2 ) ;  %  Compute  Next  Search  Point 

end; 

end; 

a  =  (max-1- ( IvuL-1 ) /2 ) /N;  %  Convert  to  Shift  Value  (Floatingpoint) 
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C.l. Log-Likelihood  Computation  Module 


LIBRARY  ieee  ; 

USE  ieee.std  logic  1164. all  ; 
USE  IEEE.std  logic  arith.all; 
USE  IEEE.std  logic  signed. all; 

ENTITY  LL  IS 
GENERIC 
( 


N 

INTEGER 

= 

12  8; 

-- 

DvL 

INTEGER 

= 

32; 

-- 

IvL 

INTEGER 

= 

32; 

-- 

DvuL 

INTEGER 

= 

3969; 

-- 

IvuL 

INTEGER 

- 

3969; 

-- 

iSg 

INTEGER 

= 

15.5; 

-- 

iSs 

INTEGER 

- 

1024; 

-- 

iSe 

INTEGER 

= 

2944; 

-- 

dSg 

INTEGER 

= 

LO 

LO 

i — 1 

-- 

dSs 

INTEGER 

- 

8; 

-- 

dSe 

INTEGER 

= 

23; 

-- 

sSs 

INTEGER 

= 

1024; 

-- 

sSe 

INTEGER 

= 

2944; 

-- 

BitDepth 

INTEGER 

= 

32 

-- 

-  Standard  Includes 


Number  of  Points  Interpolated 

Length  of  Data 

Length  of  Image 

Length  of  Data  (in  UpSampled  Terms) 

Length  of  Image  (in  UpSampled  Terms) 

Pixel  Range  of  Half  of  Image 

Start  of  Window  for  Image  (in  UpSampled  Terms) 
End  of  Window  for  Image  (in  UpSampled  Terms) 
Pixel  Range  of  Half  of  Data 
Start  of  Window  for  Data 
End  of  Window  for  Data 

First  Possible  Search  (in  UpSampled  Terms) 

Last  Possible  Search  (in  UpSampled  Terms) 
Number  of  Bits  for  Address  and  Data 


PORT 

( 


Clock,  Enable,  Reset:  IN  STD  LOGIC; 

--  Standard  Signals 

Shift 

IN 

INTEGER; 

--  Shift  Value 

Data 

IN 

INTEGER; 

--  Data  Value 

Image 

IN 

INTEGER; 

--  Image  Value 

1 Image 

IN 

INTEGER; 

--  Log  of  Image  Value 

DataAddr 

OUT 

STD  LOGIC 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Address  for  Data 

ImageAddr 

OUT 

STD  LOGIC" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Address  for  Image 

1 ImageAddr 

OUT 

STD  LOGIC" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Address  for  Log  of  Image 

Result 

OUT 

INTEGER; 

--  Log-Likelihood  Result 

Valid 

OUT 

BOOLEAN 

--  =  TRUE  When  Finished 

=  ImageAddr) 


Appendix  C:  VHDL  Version  of  Simplified  ML  Algorithm 


END  LL; 


ARCHITECTURE  extmem  OF  LL  IS 
BEGIN 

loglikelihood : 

PROCESS  (Clock,  Reset) 

VARIABLE  Idx  :  INTEGER  :=  0; 
VARIABLE  Ans  :  INTEGER  :=  0; 

BEGIN 

IF  Reset  =  ' 1 '  THEN 
Valid  <=  FALSE; 

Result  <=  -1; 

Ans  : =  -1 ; 

DataAddr  <=  (OTHERS  =>  ' 1 ' ) ; 
ImageAddr  <=  (OTHERS  =>  '  1 '  )  ; 
UmageAddr  <=  (OTHERS  =>  '  1 '  )  ; 

Idx  :  =  0 ; 

ELSIF  Rising_Edge (Clock)  THEN 
IF  Enable  =  ' 1 '  Then 

IF  (Idx  +  dSs )  >  dSe  +  1  THEN 

Valid  <=  TRUE; 

Result  <=  Ans; 

Ans  :=  Ans; 

DataAddr  <=  (OTHERS  =>  ' 1 ' ) ; 
ImageAddr  <=  (OTHERS  =>  '  1 '  )  ; 
UmageAddr  <=  (OTHERS  =>  '  1 '  ) 
Idx  :=  Idx; 

ELSE 

IF  Idx  =  0  THEN 

Valid  <=  FALSE; 

Result  <=  -1; 

Ans  :=  0; 


Computes  one  Log-Likelihood 

Loop  Variable 
Accumulates  Answer 

--  Reset  Calculator 


Displays  All 

II 

Displays  All 


1 ' s  to  Indicate  Incorrect  Answer 


1's  to  Indicate  Incorrect  Address 


Resets  Counter  for  Use 
Run  For  Loop  on  Clock 
--  Compute  Only  if  Enabled 

--  End  State,  Holds  Current  Result 

--  Outputs  Answer 

--  Displays  All  1's  to  Indicate  Incorrect  Address 


--  Stays  in  "Hold"  state 
--  Continue  Computing  Log-Likelihood 
--  If  the  Computation  Just  Started 


--  Displays  All  1's  to  Indicate  Incorrect  Answer 
--  Resets  Ans  for  Accumulation  Use 
DataAddr  <=  CONV_STD_LOGIC_VECTOR (dSs  +  idx, BitDepth) ; 

ImageAddr  <=  CONV_STD_LOGIC_VECTOR ( iSs  +  idx*N  +  (IvuL-l)/2  -  (Shift  -  1) , BitDepth) ; 
UmageAddr  <=  CONV_STD_LOGIC_VECTOR ( iSs  +  idx*N  +  (IvuL-l)/2  -  (Shift  -  1) , BitDepth) ; 

Idx  :=  Idx  +1;  --  Increment  for  Next  Clock  Cycle 

ELSE  --  If  in  Middle  of  Computation 

Ans  :=  Ans  +  Data  *  1 Image  -  Image;  --  *****  Compute  and  Add  Current  Value  to  Ans  ***** 
IF  (Idx  +  dSs)  >  dSe  THEN  --  If  the  Current  Value  is  the  End 


Valid  <=  TRUE; 

Result  <=  Ans;  --  Early  Output  of  Answer  (Saves  1  Clock  Cycle) 

DataAddr  <=  (OTHERS  =>  '1');  --  Displays  All  l's  to  Indicate  Incorrect  Address 

ImageAddr  <=  (OTHERS  =>  '  1  '  )  ;  —  " 

UmageAddr  <=  (OTHERS  =>'!');  —  " 

ELSE  --  Continues  Computation 


Valid  <=  FALSE; 

Result  <=  -1;  --  Displays  All  l's  to  Indicate  Incorrect  Answer 

DataAddr  <=  CONV_STDJLOGIC_VECTOR (dSs  +  idx, BitDepth) ; 

ImageAddr  <=  CONV_STDJLOGIC_VECTOR ( iSs  +  idx*N  +  (IvuL-l)/2  -  (Shift  -  1 ), BitDepth) ; 
UmageAddr  <=  CONV_STD_LOGIC_VECTOR ( iSs  +  idx*N  +  (IvuL-l)/2  -  (Shift  -  1), BitDepth) 


END  IF; 

I  dx  :  =  I  dx  +  1  ; 

END  IF; 

END  IF; 

ELSE 

Valid  <=  FALSE; 

Result  <=  -1; 

Ans  :=  -1; 

DataAddr  <=  (OTHERS  =>  ' 1 '  )  ; 

ImageAddr  <=  (OTHERS  =>  ' 1 ' )  ; 

UmageAddr  <=  (OTHERS  =>  '  1 '  )  ; 
Idx  :=  0; 


--  Increment  for  Next  Clock  Cycle 

--  Reset  Calculator  for  Next  Computation 

--  Displays  All  l's  to  Indicate  Incorrect  Answer 

_ fl 

--  Displays  All  l's  to  Indicate  Incorrect  Address 

_ If 

_ II 

--  Resets  Counter  for  ReUse 


END  IF; 

END  IF; 

END  PROCESS  loglikelihood; 
END  extmem; 


C.2.Main  Search  Module 


LIBRARY  ieee  ;  --  Standard  Includes 

USE  ieee.std  logic  1164. all  ; 

USE  IEEE.std  logic  arith.all; 

USE  IEEE.std  logic  signed. all; 


ENTITY  MLIW  IS 
GENERIC 


( 


N 

INTEGER 

=  12  8; 

DvL 

INTEGER 

=  32; 

IvL 

INTEGER 

=  32; 

DvuL 

INTEGER 

=  3969; 

IvuL 

INTEGER 

=  3969; 

iSg 

INTEGER 

=  15.5; 

iSs 

INTEGER 

=  1024; 

iSe 

INTEGER 

=  2944; 

dSg 

INTEGER 

=  15.5; 

dSs 

INTEGER 

=  8; 

dSe 

INTEGER 

=  23; 

sSs 

INTEGER 

=  1024; 

sSe 

INTEGER 

=  2944; 

Step 

INTEGER 

=  960; 

Grid 

BOOLEAN 

=  FALSE; 

BitDepth 

INTEGER 

=  32 

)  ; 

PORT 


--  Number  of  Points  Interpolated 

--  Length  of  Data 

--  Length  of  Image 

--  Length  of  Data  (in  UpSampled  Terms) 

--  Length  of  Image  (in  UpSampled  Terms) 

--  Pixel  Range  of  Half  of  Image 

--  Start  of  Window  for  Image  (in  UpSampled  Terms) 
--  End  of  Window  for  Image  (in  UpSampled  Terms) 

--  Pixel  Range  of  Half  of  Data 
--  Start  of  Window  for  Data 
--  End  of  Window  for  Data 

--  First  Possible  Search  (in  UpSampled  Terms) 

--  Last  Possible  Search  (in  UpSampled  Terms) 

--  Step  Size  for  Grid  Search 
--  Determines  if  Grid  Search  Happens 
--  Number  of  Bits  for  Address  and  Data 


Clock,  Enable,  Reset: 

IN  STD  LOGIC; 

--  Standard  Signals 

Data 

IN 

STD 

LOGIC 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Data  Value 

Image 

IN 

STD 

"logic" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Image  Value  (Multiplied  by  N  or  more) 

1 Image 

IN 

STD 

"logic" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Log  of  Image  Value  (Multiplied  by 

N+ ) 

DataAddr 

OUT 

STD 

"logic" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Address  for  Data 

ImageAddr 

OUT 

STD 

"logic" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Address  for  Image 

UmageAddr 

OUT 

STD 

"logic" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Address  for  Log  of  Image  (=  ImageAddr) 

Result 

OUT 

STD 

"logic" 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

--  Result:  Shift  =  (Result-1- ( IvuL-1 

/ 2 )  /N 

Valid 

OUT 

STD 

"logic" 

--  =  TRUE  When  Finished 

END  MLIW; 


ARCHITECTURE  extmem  OF  MLIW  IS 
COMPONENT  LL 
GENERIC 


N 

INTEGER 

=  12  8; 

DvL 

INTEGER 

=  32; 

IvL 

INTEGER 

=  32; 

DvuL 

INTEGER 

=  3969; 

IvuL 

INTEGER 

=  3969; 

iSg 

INTEGER 

=  15.5; 

iSs 

INTEGER 

=  1024; 

iSe 

INTEGER 

=  2944; 

dSg 

INTEGER 

=  15.5; 

dSs 

INTEGER 

=  8; 

dSe 

INTEGER 

=  23; 

sSs 

INTEGER 

=  1024; 

sSe 

INTEGER 

=  2944; 

BitDepth 

INTEGER 

=  32 

)  ; 

PORT 


Clock,  Enable,  Reset:  IN  STD  LOGIC; 

Shift 

IN 

INTEGER; 

Data 

IN 

INTEGER; 

Image 

IN 

INTEGER; 

1 Image 

IN 

INTEGER; 

DataAddr 

OUT 

STD  LOGIC 

VECTOR (BitDepth-1 

DOWNTO 

0)  ; 

ImageAddr 

OUT 

STD  LOGIC" 

VECTOR (BitDepth-1 

DOWNTO 

0) 

UmageAddr 

OUT 

STD  LOGIC" 

VECTOR (BitDepth-1 

DOWNTO 

0) ; 

Result 

OUT 

INTEGER; 

Valid 

OUT 

BOOLEAN 

)  ; 

END  COMPONENT; 


TYPE  STATE_TYPE  IS  (Start,  PreCompute,  Compute,  Hold); 
TYPE  INT_ARRAY_TYPE  IS  ARRAY  (3  DOWNTO  0)  OF  INTEGER; 
SIGNAL  State  :  STATE_TYPE ; 

SIGNAL  Answer  :  INTEGER; 


Extensible  State  Type 
INTEGER  Array  Type 
State  Machine  Variable 

Placeholder  Result  (for  Moore  Machine) 


K> 

O 


SIGNAL 

SIGNAL 

SIGNAL 

SIGNAL 

SIGNAL 

SIGNAL 


Complete 

lEnable 

IReset 

lldx 

IResult 

lValid 


BOOLEAN ; 

STD_LOGIC; 

STD_LOGIC; 

INTEGER; 

INTEGER; 

BOOLEAN ; 


BEGIN 

LL1  :  LL 
GENERIC  MAP 
( 


N 

=> 

N, 

DvL 

=> 

DvL, 

IvL 

=> 

IvL, 

DvuL 

=> 

DvuL, 

IvuL 

=> 

IvuL, 

iSg 

=> 

iSg, 

iSs 

=> 

iSs , 

iSe 

=> 

iSe, 

dSg 

=> 

dSg, 

dSs 

=> 

dSs , 

dSe 

=> 

dSe, 

sSs 

=> 

sSs , 

sSe 

=> 

sSe, 

BitDepth 

=> 

BitDepth 

) 

PORT  MAP 

( 

Clock  =>  Clock,  Enable  =>  lEnable,  Reset  =>  IReset, 
Shift  =>  lldx. 

Data  =>  CONV_INTEGER (Data)  , 

Image  =>  CONV_INTEGER ( Image) , 

1 Image  =>  CONV_INTEGER (1 Image ) , 

DataAddr  =>  DataAddr, 

ImageAddr  =>  ImageAddr, 

UmageAddr  =>  1  ImageAddr, 

Valid  =>  lValid, 

Result  =>  IResult 


Internal  State  Transition  Variable 
Enable  for  Log-Likelihood  Calculator 
Reset  for  " 

Shift  for  " 

Result  from  " 

Validity  from  " 


Moore  State  Machine  to  Control  Calcuation 


statemachine : 

PROCESS  (Clock,  Reset) 

BEGIN 

IF  Reset  =  ' 1 '  THEN 

State  <=  Start; 

ELSIF  Rising_Edge (Clock)  THEN 
CASE  State  IS 

WHEN  Start  => 

IF  Enable  =  ' 1 '  THEN 
State  <=  PreCompute; 
END  IF; 

WHEN  PreCompute  => 

IF  Complete  THEN 
State  <=  Compute; 

END  IF; 

WHEN  Compute  => 

IF  Complete  THEN 
State  <=  Hold; 

END  IF; 

WHEN  Hold  => 

IF  Enable  =  ' 1 '  THEN 
State  <=  PreCompute; 
END  IF; 

WHEN  OTHERS  => 

State  <=  Start; 

END  CASE; 

END  IF; 

END  PROCESS  statemachine; 


--  Reset  State:  Mealy  Includes  Hold  Here 
--  Moves  Immediately  to  Start  State 


--  Start  State:  Prepares  for  Calculation 

--  Waits  until  Enable  to  Move  to  PreCompute 


--  PreCompute:  Computes  EndPoints  or  Grid  Search 
--  Waits  until  Complete  to  Move  to  Compute 


--  Compute:  Performs  Gradient-Decent  Search 
--  Waits  until  Complete  to  Move  to  Hold 


--  Hold:  Outputs  Answer,  Prepares  for  Calculation 
--  Waits  until  Enable  to  Move  to  PreCompute 


--  Should  NEVER  Happen 


validout:  --  Moore  Output  of  Valid 

WITH  State  SELECT 

Valid  <=  ' 1 '  WHEN  Hold,  --  Only  in  Hold  State 

'O'  WHEN  OTHERS; 

resultout:  --  Moore  Output  of  Result 

WITH  State  SELECT 

Result  <=  CONV_STD_LOGIC_VECTOR (Answer, BitDepth)  WHEN  Hold,  —  Only  in  Hold  State 

(OTHERS  =>  '!')  WHEN  OTHERS;  --  Displays  l's  for  Incorrect  Answer 


binarysearch : 

PROCESS  (Clock,  Reset,  lValid) 


Performs  Grid  and  Binary  Searches 


K> 


VARIABLE 

VARIABLE 

VARIABLE 

VARIABLE 

VARIABLE 

VARIABLE 

BEGIN 


Idx 

NeedR 

Needs 

LLVals 

LLIdxs 

LLMax 


INTEGER 
BOOLEAN 
BOOLEAN 
:  INT_ARRAY 
:  INT_ARRAY~ 
INTEGER; 


=  0; 

=  FALSE; 
=  FALSE; 
TYPE; 
TYPE; 


Index  for  Next  Shift 

Need  Right  Value  in  Grid  Search 

Need  Slope 

LL  Values 

LL  Indexes 

Current  Max  LL  Index 


IF  Reset  =  ' 1  ' 

THEN 

-  Resets  All 

Idx  :  = 

0; 

NeedR  := 

FALSE; 

Needs  := 

FALSE; 

LLVals  := 

(OTHERS  =>  0); 

LLIdxs  := 

(OTHERS  =>  0); 

LLMax  : = 

0; 

Complete  <= 

FALSE; 

Answer  <= 

0; 

lEnable  <= 

'O'; 

IReset  <= 

'  1 '  ; 

lldx  <= 

Idx; 

ELSIF  Rising  Edge (Clock)  THEN  - 

-  Performs  Op 

CASE  State  IS 

WHEN  Start 

=> 

--  Clear  Eve 

Idx 

:  =  s  S  s  ; 

NeedR 

:=  FALSE; 

Needs 

:=  FALSE; 

LLVals 

:=  (OTHERS  => 

0)  ; 

LLIdxs 

:=  (OTHERS  => 

0)  ; 

LLMax 

:=  0; 

Complete 

<=  FALSE; 

Answer 

<=  0; 

lEnable 

<=  '  0  '  ; 

IReset 

<=  '  0 ' ; 

lldx 

<=  Idx; 

WHEN  PreCompute  => 

—  Computes 

IF  NOT  Complete  THEN 

--  Wait  fo 

IF  NOT 

lValid  THEN 

--  Wait 

Idx 

: =  Idx; 

r  Loop  to  Complete 
for  LL  Calc  to  Complete 


NeedR 

Needs 

LLVals 


NeedR; 
Needs ; 

=  LLVals; 


K> 


LLIdxs  :=  LLIdxs; 
LLMax  :=  LLMax; 
Complete  <=  FALSE; 
lEnable  <=  '  1 '  ; 

ELSIF  lEnable  =  ' 1 '  THEN 
IF  NOT  NeedS  THEN 
IF  Idx  =  sSs  THEN 
IF  NOT  Grid  THEN 
Idx  :=  sSe; 
ELSE 


--  End  of  1  LL  computation 
--  Don't  Need  the  Slope 

--  Just  store  if  it's  first 
--  Select  Next  Point 


Idx  :=  Idx  +  Step; 

END  IF; 

NeedR  :=  TRUE; 

NeedS  :=  FALSE; 

LLVals  :=  (OTHERS  =>  IResult); 

LLIdxs  :=  (OTHERS  =>  sSs); 

LLMax  : =  0 ; 

Complete  <=  FALSE; 

lEnable  <=  'O';  --  Clear  the  LL  Calculator 


ELSE  —  Not  First  Point 

IF  LLVals (LLMax)  <  IResult  THEN  --  Result  is  Bigger 
LLVals (3  DOWNTO  0)  :=  (0  =>  LLVals (3),  1  => 

2  =>  LLVals (2),  3  =>  IResult); 

LLIdxs (3  DOWNTO  0)  :=  (0  =>  LLIdxs (3),  1  => 

2  =>  LLIdxs (2),  3  =>  Idx); 


LLMax  : =  1 ; 

IF  Idx  +  Step  <=  sSe  THEN 
NeedR  :=  TRUE; 

ELSE 

NeedR  :=  FALSE; 


END  IF; 

ELSE  —  Result 

IF  NeedR  THEN 

IF  LLMax  =  0  THEN 

LLVals (3  DOWNTO  0)  :=  (0  => 

2  =>  LLVals (2) , 
LLIdxs (3  DOWNTO  0)  :=  (0  => 

2  =>  LLIdxs (2 ) , 


is  Smaller 


LLVals (0) ,  1  => 

3  =>  IResult) ; 
LLIdxs (0) ,  1  => 

3  =>  Idx) ; 


ELSE 

LLVals (3  DOWNTO  0)  :=  (0  =>  LLVals (0),  1  => 


IResult, 

Idx, 


IResult, 

Idx, 


LLVals (1) , 


K> 

4^ 


2  =>  IResult, 
LLIdxs ( 3  DOWNTO  0)  :=  (0 

2  =>  Idx, 

END  IF; 

ELSE 

LLVals ( 3  DOWNTO  0)  :=  (0 

2  =>  LLVals (2) 
LLIdxs (3  DOWNTO  0)  :=  (0 

2  =>  LLIdxs (2 ) 

END  IF; 

NeedR  :=  FALSE; 

LLMax  :=  LLMax; 

END  IF; 

IF  NOT  Grid  THEN 
Needs 
Complete 
ELSE 

IF  Idx  >=  sSe  THEN 
IF  LLIdxs (LLMax) 

Idx  :  = 

ELSE 

Idx  :  = 

END  IF; 

Needs 
ELSE 
Idx 
Needs 
END  IF; 

Complete 
END  IF; 
lEnable 


3  =>  IResult) ; 

=>  LLIdxs (0),  1  =>  LLIdxs (1) 

3  =>  Idx) ; 


=>  LLVals (0),  1  =>  LLVals (1) 

,  3  =>  IResult) ; 

=>  LLIdxs (0),  1  =>  LLIdxs (1) 

,  3  =>  Idx) ; 


--  Complete  the  EndPoints 
FALSE; 

<=  TRUE ; 

—  Continue  with  Grid  Search 
--  Setup  for  Slope  Calc 
sSe  THEN 
LLIdxs (LLMax) -1; 

LLIdxs (LLMax) +1; 

TRUE; 

--  Compute  Next  Point 
Idx  +  Step; 

FALSE; 

<=  FALSE; 

<=  'O'; 


END  IF; 

ELSE  --  Decide  Window  on  Slope 

NeedR  :=  FALSE; 

NeedS  :=  FALSE; 

IF  LLIdxs (LLMax)  =  sSs  THEN  —  Start  of  Window 
IF  LLVals (LLMax)  <  IResult  THEN 
Idx  :=  Idx; 

LLVals (1  DOWNTO  0)  :=  (0  =>  IResult,  1  =>  LLVals (1)); 

LLIdxs (1  DOWNTO  0)  :=  (0  =>  LLIdxs (LLMax) +1 ,  1  =>  LLIdxs (1)) 


=  0; 


K> 

U/i 


DOWNTO 

DOWNTO 


LLMax 
ELSE 
Idx 

LLVals (1 
LLIdxs ( 1 
LLMax  : = 

END  IF; 

ELSIF  LLIdxs (LLMax) 
IF  LLVals (LLMax)  < 
Idx  :  = 

LLVals (1 
LLIdxs ( 1 
LLMax 
ELSE 
Idx 

LLVals (1 
LLIdxs ( 1 
LLMax 
END  IF; 

ELSE 

Idx  :  = 

IF  LLVals (LLMax)  < 


=  -1; 


0)  :=  (0  =>  LLVals (0),  1  =>  LLVals (1)); 
0)  :=  (0  =>  LLIdxs (0),  1  =>  LLIdxs (1)); 
LLMax; 

=  sSe  THEN  --  End  of  Window 
IResult  THEN 
Idx; 


DOWNTO 

DOWNTO 


DOWNTO 

DOWNTO 


0)  :  = 

0)  :  = 

1; 

-1; 

0)  :  = 

0)  :  = 

LLMax; 


(0  =>  LLVals (0),  1  =>  IResult); 

(0  =>  LLIdxs (0),  1  =>  LLIdxs (LLMax) -1); 


(0  =>  LLVals (0),  1  =>  LLVals (1)); 
(0  =>  LLIdxs (0),  1  =>  LLIdxs (1)); 


--  Middle  of  Window 

Idx; 

IResult  THEN 


LLVals (1 

DOWNTO 

0) 

:=  (0 

=> 

IResult, 

1 

=>  LLVals (2) ) ; 

LLIdxs ( 1 

DOWNTO 

0) 

:=  (0 

=> 

LLIdxs (LLMax) +1, 

1  =>  LLIdxs (2 

LLMax 

:  = 

0; 

ELSE 

LLVals (1 

DOWNTO 

0) 

:=  (0 

=> 

LLVals (0),  1 

=> 

LLVals (1) ) ; 

LLIdxs ( 1 

DOWNTO 

0) 

:=  (0 

=> 

LLIdxs ( 0 ) ,  1 

=> 

LLIdxs ( 1 ) ) ; 

LLMax 

:  = 

1; 

END  IF; 
END  IF; 


LLVals (3 

DOWNTO  2) 

LLIdxs ( 3 

DOWNTO  2) 

Complete 

<=  TRUE 

lEnable 

<=  '  0  '  ; 

END  IF; 

LSE 

Idx  :  = 

Idx; 

NeedR  := 

NeedR; 

Needs  : = 

Needs ; 

:=  LLVals (3  DOWNTO  2); 
:=  LLIdxs (3  DOWNTO  2); 


--  Wait  for  LL  Calc  to  Clear 


K> 

0\ 


LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 

LLMax  :=  LLMax; 

Complete  <=  FALSE; 
lEnable  <=  lEnable; 

END  IF; 

ELSE  --  Leaving 

IF  Idx  =  -1  THEN 
Idx  :=  Idx; 


ELSE 

Idx  :=  (LLIdxs (0) ILLIdxs (1) ) /2; 

END  IF; 

NeedR  :=  FALSE; 

NeedS  :=  FALSE; 

LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 

LLMax  :=  LLMax; 

Complete  <=  FALSE; 

lEnable  <=  'O'; 


END  IF; 

Answer  <=  0; 

IReset  <=  'O'; 

lldx  <=  Idx; 


State 


Setup  for  Compute 


WHEN  Compute  => 

IF  NOT  Complete  THEN 
IF  NOT  lValid  THEN 
Idx  :=  Idx; 

Needs  :=  NeedS; 

LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 


LLMax  :=  LLMax; 


IF  Idx  =  -1  THEN 

Complete  <=  TRUE; 
lEnable  <=  'O'; 

ELSE 


Complete 
lEnable 
END  IF; 
ELSIF  lEnable 


FALSE; 

'  1 ' ; 

' 1 '  THEN 


Performs  Full  Search 
-  Wait  for  Loop  to  Complete 
--  Wait  for  LL  Calc  to  Complete 


--  End  of  1  LL  computation 


<= 

<= 


K> 
■ — 1 


IF  NOT  Needs  THEN  --  Do  Pre-Slope  Computations 

IF  LLVals (LLMax)  <  IResult  THEN  --  Result  is  Bigger 
LLVals(2)  :=  IResult; 

LLIdxs(2)  :=  Idx; 


IF  Idx  +  1 

<  LLIdxs (1)  THEN 

-■ 

-  Setup  for  Slope  Calc 

Needs 

:  = 

TRUE  ; 

LLVals (1 

DOWNTO 

0)  :  = 

(0 

=> 

LLVals (0),  1  =>  LLVals (1) 

LLIdxs ( 1 

DOWNTO 

0)  :  = 

(0 

=> 

LLIdxs (0),  1  =>  LLIdxs (1) 

LLMax 

:  = 

LLMax; 

Idx 

:  = 

I dx  +  1 ; 

ELSE 

-■ 

-  To  Left,  Shutter  Right 

Needs 

:  = 

FALSE; 

LLVals (1 

DOWNTO 

0)  :  = 

(0 

=> 

LLVals (0),  1  =>  IResult); 

LLIdxs ( 1 

DOWNTO 

0)  :  = 

(0 

=> 

LLIdxs (0),  1  =>  Idx); 

LLMax 

:  = 

1; 

Idx 

:  = 

(LLIdxs (0) 

+ I dx ) / 2 ; 

END  IF; 

LSE 

—  ] 

Result  is  Smaller 

Needs 

:  = 

FALSE; 

LLVals (2) 

:=  LLVals (2); 

LLIdxs (2 ) 

: =  LLIdxs (2 ) ; 

IF  LLMax  = 

0  THEN 

-■ 

-  To  Left,  Shutter  Right 

LLVals (1 

DOWNTO 

0)  :  = 

(0 

=> 

LLVals (0),  1  =>  IResult); 

LLIdxs ( 1 

DOWNTO 

0)  :  = 

(0 

=> 

LLIdxs (0),  1  =>  Idx); 

Idx 

:  = 

(LLIdxs (0) 

+Idx) / 2 ; 

ELSE 

-■ 

-  To  Right,  Shutter  Left 

LLVals (1 

DOWNTO 

0)  :  = 

(0 

=> 

IResult,  1  =>  LLVals (1)) 

LLIdxs ( 1 

DOWNTO 

0)  :  = 

(0 

=> 

Idx,  1  =>  LLIdxs (1)); 

Idx 

:  = 

( Idx+LLIdxs ( 1 

)  )  /2; 

END  IF; 

LLMax 

;  = 

LLMax; 

END  IF; 

lEnable  <=  'O';  --  Clear  LL  Calculator 

ELSE  --  Do  Post-Slope  Computations 

NeedS  :=  FALSE; 

IF  LLVals (2)  <  IResult  THEN  —  To  Right,  Shutter  Left 
LLVals (1  DOWNTO  0)  :=  (0  =>  IResult,  1  =>  LLVals (1)); 

LLIdxs ( 1  DOWNTO  0)  :=  (0  =>  Idx,  1  =>  LLIdxs (1)); 

LLMax  : =  0 ; 

Idx  :=  ( Idx+LLIdxs ( 1 ) ) /2 ; 


ELSE  --  To  Left,  Shutter  Right 

LLVals ( 1  DOWNTO  0)  :=  (0  =>  LLVals(O),  1  =>  LLVals (2)); 
LLIdxs ( 1  DOWNTO  0)  :=  (0  =>  LLIdxs(O),  1  =>  LLIdxs (2)); 
LLMax  : =  1 ; 

Idx  :=  (LLIdxs (0) +LLIdxs (2) )/2; 


END  IF; 

LLVals (2)  :=  LLVals (2); 

LLIdxs (2)  :=  LLIdxs (2); 

END  IF; 

LLVals (3)  :=  LLVals (3); 

LLIdxs (3)  :=  LLIdxs (3); 

IF  Idx  >  LLIdxs (0)  AND  Idx  < 


Complete 

<= 

FALSE; 

ELSE 

Complete 

<= 

TRUE  ; 

END  IF; 

lEnable 

<= 

'O'; 

ELSE 

Idx  :  = 

Idx; 

Needs  := 

Needs ; 

LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 

LLMax  :=  LLMax; 
Complete  <=  FALSE; 

lEnable  <=  lEnable; 


LLIdxs (1)  THEN  --  Continue  Search 
--  End  of  Search 


Wait  for  LL  Calculator  to  Clear 


END  IF; 

ELSE 

Idx  :=  Idx; 

NeedS  :=  FALSE; 

LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 

LLMax  :=  LLMax; 

Complete  <=  FALSE; 

lEnable  <=  'O'; 


Leaving  State  -  Setup  for  Hold 


END  IF; 

NeedR  :=  FALSE; 

Answer  <=  0; 

IReset  <=  'O'; 

lldx  <=  Idx; 


'vO 


WHEN  Hold  =>  --  Hold  Output  Answer 

I  dx  :  =  0 ; 

NeedR  :=  FALSE; 

NeedS  :=  FALSE; 

LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 

LLMax  :=  LLMax; 

Complete  <=  FALSE; 

Answer  <=  LLIdxs (LLMax) ; 
lEnable  <=  ' 0  '  ; 

IReset  <=  'O'; 

lldx  <=  Idx; 


WHEN  OTHERS  =>  —  Should  Never  Happen 

Idx  : =  0 ; 

NeedR  :=  FALSE; 

NeedS  :=  FALSE; 

LLVals  :=  LLVals; 

LLIdxs  :=  LLIdxs; 

LLMax  :=  LLMax; 

Complete  <=  Complete; 

Answer  <=  Answer; 

lEnable  <=  'O'; 

IReset  <=  'O'; 

lldx  <=  Idx; 

END  CASE; 

END  IF; 

END  PROCESS  binarysearch; 


END  extmem; 


Appendix  D:  Complete  Parameterization  of  Bias  and  Noise  Statistics 


D.l.LGS  Image  Size  Comparison 


Parameters:  C  =  300,  O;  =  2,  Bg  =  0,  lxNyquist 


Minimum  Sizes 


Recommended  Sizes 
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MSE  vL&S,  1/2  Wave  of  Tilt  MSE  v  L  &  Bg,  1/2  Wave  of  Tilt  MSE  v  L  &  C,  1/2  Wave  of  Tilt 


D.2.LGS  swat  Characterization 


Parameters:  Image  Size  (L),  Intensity  (C)  =  300,  a;  =  2,  Background  (Bg)  =  0,  &  Nyquist 
Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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MSE  vL&S,  1/2  Wave  of  Tilt  MSE  v  L  &  Bg,  1/2  Wave  of  Tilt  MSE  v  L  &  C,  1/2  Wave  of  Tilt 


D.3.LGS  mliwi  Characterization 


Parameters:  Image  Size  (L),  Intensity  (C)  =  300,  a;  =  2,  Background  (Bg)  =  0,  &  Nyquist 
Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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MSE  vL&S,  1/2  Wave  of  Tilt  MSE  v  L  &  Bg,  1/2  Wave  of  Tilt  MSE  v  L  &  C,  1/2  Wave  of  Tilt 


D.4.LGS  mliwc  Characterization, 


Parameters:  Image  Size  (L),  Intensity  (C)  =  300,  a;  =  2,  Background  (Bg)  =  0,  &  Nyquist 
Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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D.5.  Tracking  Hubble  Image  Size  Comparison. 

Parameters:  C=8000,  a,~2,  Bg=0,  lxNyquist 


Minimum  Sizes 


Recommended  Sizes 
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Shift  Value  {pixels}  Shift  Value  {pixels}  Shift  Value  {pixels} 


D.6.  Tracking  Hubble  mliwi  Characterization  (y  dim). 

Parameters:  Image  Size  (L),  Intensity  (C)  =  8000,  o,=  2,  Background  (Bg)  =  0,  & 
Nyquist  Sampling  (S)  =  l,  Unless  Otherwise  Noted. 
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L  {pixels}  1  S  {unitless}  L  {pixels}  0  Bg  {photons}  L  {pixels}  C  {photons} 


D.7.  Tracking  Hubble  mliwc  Characterization  (y  dim). 

Parameters:  Image  Size  (L),  Intensity  (C)  =  300,  a;  =  2,  Background  (Bg)  =  0,  &  Nyquist 
Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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Absolute  Bias  v  Shift  VAR  v  Shift  WISE  v  Shift 


D.8.WFS  Hubble  Image  Size  Comparison 


Parameters:  C=300,  a,~2,  Bg=0,  lxNyquist 


Minimum  Sizes 


Recommended  Sizes 
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Shift  Value  {pixels}  Shift  Value  {pixels}  Shift  Value  {pixels} 


D.9.WFS  Hubble  swat  Characterization 

Parameters:  Image  Size  (L),  Intensity  (C)  =  8000,  a;  =  2,  Background  (Bg)  =  0,  & 
Nyquist  Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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D.10.  WFS  Hubble  mliwi  Characterization 

Parameters:  Image  Size  (L),  Intensity  (C)  =  8000,  a;  =  2,  Background  (Bg)  =  0,  & 
Nyquist  Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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L  {pixels}  1  S  {unitless}  L  {pixels}  0  Bg  {photons}  L  {pixels}  C  {photons} 


MSE  vL&S,  1/2  Wave  of  Tilt  MSE  v  L  &  Bg,  1/2  Wave  of  Tilt  MSE  v  L  &  C,  1/2  Wave  of  Tilt 


D.ll.  WFS  Hubble  mliwc  Characterization 


Parameters:  Image  Size  (L),  Intensity  (C)  =  8000,  a;  =  2,  Background  (Bg)  =  0,  & 
Nyquist  Sampling  (S)  =  1,  Unless  Otherwise  Noted. 
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